NASA/TM-2014-218243

Modeling and Simulation of Upset-Inducing 

Disturbances for Digital Systems in an 

Electromagnetic Reverberation Chamber

Wilfredo Torres-Pomales

Langley Research Center, Hampton, Virginia

April 2014
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

The NASA STI program operates under the 

auspices of the Agency Chief Information Officer. 

It collects, organizes, provides for archiving, and 

disseminates NASA’s STI. The NASA STI 

program provides access to the NASA Aeronautics 

and Space Database and its public interface, the 

NASA Technical Report Server, thus providing one 

of the largest collections of aeronautical and space 

science STI in the world. Results are published in 

both non-NASA channels and by NASA in the 

NASA STI Report Series, which includes the 

following report types: 

 TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

Programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but 

having less stringent limitations on manuscript 

length and extent of graphic presentations. 

 TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

 CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

 CONFERENCE PUBLICATION. 

Collected papers from scientific and 

technical conferences, symposia, seminars, 

or other meetings sponsored or cosponsored by NASA. 

 SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, 

often concerned with subjects having 

substantial public interest. 

 TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission. 

Specialized services also include organizing 

and publishing research results, distributing 

specialized research announcements and feeds, 

providing information desk and personal search 

support, and enabling data exchange services. 

For more information about the NASA STI 

program, see the following: 

 Access the NASA STI program home page 

at http://www.sti.nasa.gov

 E-mail your question to help@sti.nasa.gov

 Fax your question to the NASA STI 

Information Desk at 443-757-5803 

 Phone the NASA STI Information Desk at 

443-757-5802 

 Write to: 

 STI Information Desk 

 NASA Center for AeroSpace Information 

 7115 Standard Drive 

 Hanover, MD 21076-1320 
NASA/TM-2014-218243

Modeling and Simulation of Upset-Inducing 

Disturbances for Digital Systems in an 

Electromagnetic Reverberation Chamber

Wilfredo Torres-Pomales

Langley Research Center, Hampton, Virginia

National Aeronautics and 

Space Administration

Langley Research Center 

Hampton, Virginia 23681-2199

April 2014
Available from: 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 

443-757-5802 

Acknowledgment 

The author is grateful for the comments and suggestions by Dr. Van E. Brewer of Old 

Dominion University, and Kenneth W. Eure, Truoung X. Nguyen, Laura J. Smith and 

Allan L. White of NASA Langley Research Center. 
Abstract

This report describes a modeling and simulation approach 

for disturbance patterns representative of the environment 

experienced by a digital system in an electromagnetic 

reverberation chamber. The disturbance is modeled by a multivariate statistical distribution based on empirical observations. 

Extended versions of the Rejection Samping and Inverse 

Transform Sampling techniques are developed to generate multivariate random samples of the disturbance. The results show 

that Inverse Transform Sampling returns samples with higher 

fidelity relative to the empirical distribution. This work is part of 

an ongoing effort to develop a resilience assessment 

methodology for complex safety-critical distributed systems.

iii
Table of Contents

1. Introduction ................................................................................................................................................ 1

2. Modeling Approach.................................................................................................................................... 7

3. Simulation Approach.................................................................................................................................. 8

4. Results .......................................................................................................................................................12

4.1. Empirical Probability Distribution Function ......................................................................................12

4.2. Pseudo-Random Number Generator...................................................................................................12

4.3. Overstroke Sample Generation...........................................................................................................14

4.4. Goodness-of-Fit Measures..................................................................................................................17

5. Conclusions ...............................................................................................................................................18

References .....................................................................................................................................................18

Abbreviations ................................................................................................................................................20

iv
1. Introduction

A research effort is underway to develop practical validation and verification (V&V) methods that can 

enable rigorous safety assurance for the next generation of aviation systems. T hese systems are 

characterized by highly complex, large-scale, network-based distributed architectures with softwareimplemented functionality and a dvanced computation and communication capabilities. To meet the 

safety goals, these systems must be dem onstrably robust w ith respect to system design and 

implementation errors, component degradations and failures, and partial system failures. The V&V 

challenge is compounded by strong coupling of system components in the software and the hardware, as 

well as the need to consider unexpected and possibly malicious component behavior [1]. 

To support this research effort, verification approaches are being developed for robust distributed 

algorithms that support system resource management in a fault space with a wide range of severity. A 

system architecture for safety-critical real-time applications must have the ability to mitigate the effects of 

internal component faults of varying severity [2]. A safety-critical system must have sufficient tolerance 

to accommodate more frequent uncorrelated random faults without disruption of the system services. A 

robust system must also mitigate infrequent but more severe correlated faults that can exceed the system 

design assumptions, disrupt internal coordinated operation among the system components and propagate 

effects outward to the external service interfaces. Analysis techniques will be developed for system 

designs intended to en sure continued safe operation in the p resence of c omponent misbehavior while 

simultaneously minimizing their adverse effects. These techniques should enable designs with strongly 

assured safety properties under the weakest possible (i.e., least restrictive) assumptions in terms of the 

number and types of faults a system can handle. 

In this context, we are developing analysis and simulation tools to assess and gain insight into the error 

propagation dynamics of fault-tolerant distributed architectures. In [3] we proposed an approach to assess 

system resilience, which we defined as the ability to contain and recover from fault effects. The threats

to the delivery of proper system service are internal system conditions or external environmental 

conditions (e.g., HIRF, lightning, high-energy particle radiation, power system transients, etc.) that may 

cause faults in the system. There is special interest in examining the response to functional upsets, which 

are error modes that involve no pe rmanent component damage, can s imultaneously occur in multiple 

channels of a redundant distributed system and can cause unrecoverable state error conditions [4, 5, 6].

Under the proposed assessment approach, the system is viewed from a stimulus-response (i.e., cause 

and effect) perspective (see Figure 1 and Figure 2). The threat conditions specify the stimulus space,

which is a subset of all possible system threat patterns. A disturbance is an external system stimulus that 

may cause a perturbation, defined here as an internal fault condition in the form of outages on the 

services provided by the components. Alternatively, the specification of the disturbance can be bypassed 

and the stimulus can be specified as a perturbation. The e ffects of a perturbation (i.e., errors) may 

propagate throughout the system and reach the external functional interface, thus causing an outage of the 

external system service, which we refer to as a disruption. The response space is the set of system 

disruptions resulting from the app lication of the stimulus space. A s Figure 2 suggests, for a g iven 

perturbation, the severity of the disruption is determined by the error propagation characteristics of the 

system.

1
For a quantitative system resilience analysis, the disturbance, perturbation and disruption spaces may 

be described in terms of probability distributions (PDs) of random variables whose values correspond to 

the severity of occurrences (i.e., events or instances) in these spaces. Figure 3 illustrates the use of PDs to 

describe the spaces in the stimulus-response chain. To enable the use of such distributions in analyses, 

we need severity metrics for the occurrences in each space. A simple occurrence severity model for a 

HIRF disturbance was described in [7]. Fo r perturbations and disruptions, we use the concept of 

corruption, which is defined as the amount of error in a service outage. Our proposed corruption metrics 

for perturbations and disruptions are presented in [3]. Given the probability distribution for a stimulus, 

we can perform an injection experiment (e.g., a Monte Carlo experiment) to generate a response set for 

which we can then compute the probability distribution. The design goal for a resilient system is to 

mitigate the fault space by creating an inv erse relation between the probability of occur rence and the

severity of effects.

As part of this effort to d evelop robust resilient systems, two exper iments were condu cted at th e 

Component-Level Interfaces: 

Internal Structural Services

Components

System

System-Level Interface:

Functional Service

External 

or Internal 

Threats

Stimulus Response

Perturbation:

Internal Response 

and Stimulus

Disruption:

External 

Response

Disturbance:

External 

Stimulus

Failure

Units

Error 

Propagation

System

Probability

Severity

Disturbance Space

Severity

Perturbation Space

Severity

Disruption Space

Failure

Units

System 

Structure

Probability Probability

Figure 1: Stimulus-response system model

Figure 2: Stimulus-response chain

Figure 3: Description of spaces in the stimulus-response chain using probability distributions

2
NASA Langley High Intensity Radiated Field Laboratory [8] to investigate the effects of radiated 

electromagnetic fields on digital computer systems. During the first experiment, a prot otype of an 

onboard data communication system was exposed to a wide range of electromagnetic field conditions to 

measure the effects of the radio frequency (RF) radiation on the performance of the system [9]. Figure 4

shows a typical observed system response. Every time an RF radiation burst was applied to the system, 

the same or similar response pattern was observed. This repetitive response led to speculation that the 

radiation in the chamber also had a repetitive pattern1

. Figure 5 shows a sample of the peak-normalized 

electromagnetic field strength measured by an antenna inside the chamber. The field amplitude pattern 

was indeed found to be repetitive and dependent on chamber conditions like frequency, input power and 

others. Based on this, it was conjectured that the error burst pattern was caused by the pattern of field 

excursions above the susceptibility threshold of the radiated device. We refer to these excursions of the 

field strength amplitude as overstrokes. Figure 7 shows an overstroke pattern generated by “slicing” a 

field amplitude curve at a level representing the su sceptibility threshold of a d evice2

. The overstroke 

amplitude of 0.0 corresponds to the field amplitude at the susceptibility threshold. We believe that the 

error rate and duration of functional error bursts are strongly correlated with the amplitude and duration of 

field strength excursions above the susceptibility threshold of the radiated device. This error generation 

mechanism could help explain the relation between field strength and error probability reported by Yates 

et al. [10]. 

A second experiment was conducted to collect field amplitude data suitable for a characterization of 

the electromagnetic environment in a reverberation chamber [7]. The characterization of the 

electromagnetic field was based on statistics of overstroke peak amplitude and duration for a range of 

susceptibility threshold profiles. Am ong the repor ted analyses, we c onsidered the d istribution of 

overstrokes for the susceptibility threshold profile measured in the first experiment. Figure 6 shows the 

scatter plot for the corresponding overstrokes. There are 4340 points in Figure 6. Figure 8 shows the 

histograms of the relative frequency distribution of overstroke peak amplitude and duration when the 

overall peak field strength is 20 V/m above the susceptibility threshold (Eos = 20 V/m).

By combining the results from these two experiments, we could develop a simple fault simulation 

model in which ran domly generated disturbances in the form of field overstrokes are mapped to error 

bursts on system models. This relation between the disturbances and perturbations is a critical element in 

the approach described above for the characterization of sy stem resilience with respect to a pa rticular

disturbance space [3]. 

The goal of the work presented in th is report was to dev elop a rand om sample generator for an 

overstroke population described in terms of the peak and dur ation of the overstrokes. We w ant to 

generate random samples that a re not t aken directly from the g iven population sam ple but hav e 

approximately the sam e statistical properties. We a lso want t he sample generation technique to be 

generic but easily tunable so it can be applied to any given overstroke population. To meet this goal, we 

applied a generic, non-parametric statistical model for field overstrokes and a simple simulation solution 

to generate random samples from this model. 

1 This is actually a well-known fact due to the repeated patterns of the stirrers, but we were not aware of it at the 

time of the experiment.

2 Note that the data in Figures 4, 5, and 6 are for illustrative purposes only and they are not directly related. The data 

in these figures were not taken at the same time or with the same test conditions.

3
Normalized Field Strength Trace

Frequency: 121.15 MHz

0.00

0.10

0.20

0.30

0.40

0.50

0.60

0.70

0.80

0.90

1.00

1.10

0.000.130.250.380.510.630.760.881.011.141.261.391.521.641.771.892.022.152.272.402.532.652.782.903.033.163.283.413.543.663.793.914.044.174.294.424.554.674.804.92

Stirrer Position (sec)

Normalized Field Strength

Error Rate During a Radiation Strike

Frequency: 121.15 MHz; Overstrength: 20 V/m

0.00

0.10

0.20

0.30

0.40

0.50

0.60

0.70

0.80

0.90

1.00

14.0014.1014.2014.3014.4014.5014.6014.7014.8014.9015.0015.1015.2015.3015.4015.5015.6015.7015.8015.9016.0016.1016.2016.3016.4016.5016.6016.7016.8016.9017.0017.1017.2017.3017.4017.5017.6017.7017.8017.9018.0018.1018.2018.3018.4018.5018.6018.7018.8018.9019.00

Time (sec)

Error Rate

Overstrokes Susceptibility Threshold Level = 0.85

Figure 4: Example of system response to RF radiation

Figure 5: Sample trace of the electromagnetic field strength in a reverberation chamber

4
Figure 6: Scatter plot of normalized experimental overstroke peak amplitude and duration

Stirrer-Induced Modulation Overstrokes

Frequency: 121.15 MHz; Reference Level: 0.65

0.000

0.050

0.100

0.150

0.200

0.250

0.300

0.350

0.400

200.00200.10200.20200.30200.40200.50200.60200.70200.80200.90201.00201.10201.20201.30201.40201.50201.60201.70201.80201.90202.00202.10202.20202.30202.40202.50202.60202.70202.80202.90203.00203.10203.20203.30203.40203.50203.60203.70203.80203.90204.00204.10204.20204.30204.40204.50204.60204.70204.80204.90205.00

Time (sec)

Overstroke Amplitude 0.00

0.10

0.20

0.30

0.40

0.50

0.60

0.70

0 100 200 300 400 500 600 700 800 900 1000

Peak Amplitude (normalized)

Duration (ms)

Distribution of Overstroke Peak Amplitude and Duration

Nominal Overstrength = 20.0 V/m

Figure 7: Example of overstrokes for a particular susceptibility threshold level

5
Figure 8: Histograms of experimental overstroke amplitude and duration

6
The next section describes the selected approach. This is followed by a description of the results. The 

report ends with brief concluding remarks summarizing what has been accomplished and future work.

2. Modeling Approach

A primary constraint in this study was that only the existing experimental field data as reported in [7]

was available for modeling and si mulation activities. It was possible to generate samples of different 

overstroke populations with the existing data (for example, by changing the susceptibility threshold used 

to determine the overstrokes from the exi sting field amplitude traces), but there was only one set o f 

overstrokes for any given population as defined by the param eter configuration used to g enerate the 

overstrokes from field amplitude traces. This could be an issue if we tried to assess the absolute accuracy 

of the dev eloped model. Since t his was no t a hy pothesis testing study [11], the assessment of 

representativeness of a g iven experimental overstroke sample was no t within the scope of the study.

Rather, the problem was framed such that we had to determine how well t he distribution of simulated 

overstrokes matched a given experimental sample distribution (i.e., compare experimental vs. simulated 

samples), assuming that the given sample was representative of the disturbance pattern that an actual 

system would experience in a reverberation chamber. Only the overstroke sample shown in Figure 6 was 

considered in this study.

An overstroke population is characterized in terms of two continuous random variables, the peak 

amplitude P and the duration D, and a b ivariate probability distribution. An overstroke observation 

consists of a pa ir (p, d) c orresponding to th e peak amplitude and duration of the ov erstroke. For th e 

empirical sample, P and D have ranges [0, pmax] and [0, dmax].

The first modeling approach considered was the use of em pirical marginal cumulative distribution 

functions (CDF) for the variables P and D, and an empirical copula3 to capture the dependence between P

and D in terms of a joint cumulative distribution function for the marginal CDFs [12]. Although simple 

in concept, this approach was more complex to model and simulate than the alternative approaches. 

The second alternative modeling approach was an empirical bivariate cumulative distribution function 

(eCDF), which is defined as [13]:

FP,D(p, d) = [6I I(pi d p, di d d)]/n,

for an n-point sample where the i- th point is denoted (pi, di). I(x, y) is the indi cator function, which is 

equal to 1 when x and y are true, and zero otherwise. This simple modeling approach would work well 

for our sample data. However, it was rejected because the simulation based on it was more complex than 

for our next alternative.

The selected modeling approach was an em pirical probability distribution function (ePDF). For the 

bivariate case, this is simply a three-dimensional histogram with rectangular bins on the horizontal plane 

and the vertical dimension equal to the relative frequency of experimental points in each bin. As for the 

basic two-dimensional histogram, this approach works well if the sample size is large. Given that we had 

4340 points in the experimental sample and most of them were concentrated close to the (0, 0) origin (see 

3 A copula is a multivariable probability distribution in which the marginal distributions are uniform.

7
Figure 6), the ePDF was the easiest and most effective alternative modeling approach. Additionally, as 

described in the next section, there were at least two simple techniques that could be applied to generate 

random samples based on the ePDF. 

To build a bivariate ePDF, we had to specify the number and width of the bins. Two options were 

considered for the number of bins [14]. The selected option was the square-root rule, where the number 

of bins is given by:

k = ξn.

With this choice, the growth in the size of the model (i.e., its complexity) is sub-linear relative to the 

number of sample points. 

For the bin width selection to work in combination with the square-root rule for the number of bins, 

we considered two alternatives: constant bin widths and percentile-based bin widths. For the constant bin 

width option, the bin widths for P and D, denoted hp and hd, are given by:

hp = pmax/k and hd = dmax/k.

The percentile-based bin width option is intended to give the model a higher resolution around the 

intervals with the highest concentration of points. For this bin-width alternative, the upper bound of the 

bin corresponding to the p-th percentile for the P random variable is given by the i-th value of P when 

arranged in increasing order:

i = Round(p(n+ 1)/100).

The bin widths are determined independently for the ranges of the P and D random variables.

The results presented in Section 4 show t hat the empirical probability distribution function using the 

square-root rule and constant bin widths is an adequate model for our experimental sample distribution.

Alternatively, we can use the Freedman-Diaconis rule with constant bin width given by:

h = 2 (IQR/n1/3)

where IQR is the int er-quartile range (i.e., the distance between the 25th and 75th percentiles). This bin 

width takes into consideration the actual distribution of the sample and is insensitive to outliers. This rule 

would be app lied independently to the peak P and distribution D rand om variables to d etermine the 

dimensions of the bins. The Freedman-Diaconis rule is available for future use, if needed, but it was not 

implemented in the model and simulations presented in this report.

3. Simulation Approach

Simulation is simply the dynamic generation of samples from the modeled population. This is realized 

by a mathematical transformation based on the modeled population that takes uniform random variables 

as inputs and g enerates random variables with the desired distribution. One or m ore pseudo-random 

number generators (PRNGs) are used to produce samples of uniformly distributed variables. There are a

8
number of techniques that can be used to transform uniform random samples into the desired population, 

including Rejection Sampling, Inverse Transform Sampling, Metropolis-Hasting algorithm, Gibbs 

sampling, Slice sampling, and others [15]. T he main difficulty for our simulation problem is that the 

random variables that characterize the overstroke population, P and D, are not independent, and thus, we 

need a way to accurately capture their correlation. In addition, we wanted our simulation approach to 

require only a few simple mathematical operations, be easy to implement and be computationally fast.

We considered two sampling techniques: Rejection Sampling and Inverse Transform Sampling.

Rejection Sampling (also known as A cceptance-Rejection Method) can be easily understood using 

Figure 9. Suppose we had a random variable r with range [0, 1] and PDF g iven by the red l ine. Let u1 

and u2 denote two uniform and independent (i.e., uncorrelated) random variables, both with range [0, 1]. 

If we generate uniformly distributed random samples in the u1 x u2 sp ace and we discard all samples

above the r ed line, the samples that remain are distributed according to the P DF of r. I n effect, this 

algorithm transforms u1 and u2 into r. 

To apply Rejection Sampling to the overstroke bivariate distribution, we transform three uniform and 

uncorrelated random variables U1, U2, U3 into correlated variables P and D using the empirical PDF as the 

rejection threshold. Specifically, sample point (p, d) is given by:

p = pmax·u1, and

d = dmax·u2.

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0 0.2 0.4 0.6 0.8 1

u2

u1

Figure 9: Illustration of Rejection Sampling concept

Rejected

Accepted

9
Let fP,D(p, d) denote the overstroke PDF a t point (p, d). I f u3 ! fP,D(p, d), point (p, d) is rejected; 

otherwise, (p, d) is accepted. Comparison with fP,D(p, d) ensures that P and D are properly correlated.

Application of the Rejection Sampling technique requires the evaluation of fP,D(p, d) at every point in 

the P x D space as sample points (p, d) computed from variables U1 and U2 have continuous, real-valued 

ranges. We used the bilinear interpolation technique to extend our ePDF model from a discrete domain to 

a continuous domain. We use Figure 10 to illustrate the application of bilinear interpolation to determine 

the value of Q = fP,D(p, d). In Figure 10, we are operating in a three-dimensional space and we have four 

known reference points defined by our ePDF model: (p1, d1, Q11), (p2, d1, Q21), (p1, d2, Q12), and (p2, d2,

Q22). We assume that p1 < p < p 2 and d1 < d < d 2. To determine Q, we first interpolate along the 

horizontal dimension in Figure 10 and then along the vertical dimension. We define four scaling factors:

wp1 = (p – p1)/(p2 – p1)

wp2 = (p2 – p)/(p2 – p1)

wd1 = (d – d1)/(d2 – d1)

wd2 = (d2 – d)/(d2 – d1)

We interpolate along the horizontal dimension to determine R1 and R2:

R1 = wp2Q11 + wp1Q21

R2 = wp2Q12 + wp1Q22

Next, we interpolate along the vertical dimension to determine Q:

Q = wd2R1 + wd1R2 = wp2wd2Q11 + wp1wd2Q21 + wp2wd1Q12 + wp1wd1Q22

The corresponding equations for Q when p = p2 and/or d = d2 require interpolation along at most one 

dimension and are easily derived.

x x

x x

x

x

x

x x

Q11

x

Q21

x

Q12 Q22

x

p1 p

x

p2

x

Q

p1

d1

d x

x

Q

d2

x

Q

x

R1

x

R2

x

S2

x

S1

Figure 10: Reference points of bilinear interpolation technique

10
Our second simulation technique was a form of Inverse Transform Sampling applied to a bivariate 

distribution. Inverse Transform Sampling relies on the Probability Integral Transform (PIT) [16]. Let P 

and U denote r andom variables, and FP(P) denote the marginal CDF of P. If U = FP(P), then U has a 

uniform distribution with range [0, 1]. For the Inverse Transform Sampling technique, we apply the PIT 

transform in the opposite direction: P = FP

-1(U). This is illustrated in Figure 11 where we use the random 

value u to graphically determine the corresponding value p. T o expand the inverse transform to the 

bivariate case, we use two uniform random variables U1 and U2. First, we apply the inverse transform to 

the marginal CDF of P:

p = FP

-1(u1)

Then, we apply the inverse transform to the conditional CDF of d, denoted FD|P(d|p):

d = FD|P

-1(u2|p)

The marginal FP(P) is simply a cumulative histogram of the overstroke sample set considering only the P 

values. The number of bins was set using the same square-root rule as for the ePDF model. For the 

conditional CDF of d, F D|P(d|p), we appli ed bilinear interpolations to the ePDF to det ermine the 

conditional PDF of d given p, denoted fD|P(d|p), and then computed the conditional CDF, FD|P(d|p). 

0

0.2

0.4

0.6

0.8

1

1.2

0.000

33.12366.24599.368

132.491165.614198.736231.859264.982298.105331.227364.350397.473430.595463.718496.841529.964563.086596.209629.332662.455695.577728.700

Cumulative Relative Frequency

p

x

u

p

6

8

x

4568

p

x

Figure 11: Illustration of Inverse Transform Sampling

11
4. Results

The results cover the aspects of modeling and simulation, including pseudo-random number 

generation, and comparison of R ejection Sampling and I nverse Transform Sampling. We created C 

language programs to produce the ePDF model and simulate it to generate random samples.

4.1. Empirical Probability Distribution Function

Figure 12 shows the ePDF using the square-root rule and constant bin widths. The results presented 

below show that this is an adequate model for our experimental sample distribution.

4.2. Pseudo-Random Number Generator

The simulation of the model requires pseudo-random number generators (PRNGs) whose output are

approximately uniform and uncorrelated between samples and between generators. Figure 13 shows the 

PDF for a PRNG. As expected, each of the 10 bins contains approximately one-tenth of the generated 

random numbers. Fig ure 14 shows that the cumulative distribution function of t he PRNGs g rows 

linearly, as exp ected for a uniform random variable. Fig ure 15 shows a scatter plot where the point 

coordinates are the outputs of two PRNGs: (PRNG1, PRNG2). Thi s figure shows that there is n o 

discernible pattern of correlation between the PRNGs. These results provide confidence that the PRNGs, 

as the foundation for our simulation approach, meet the bas ic requirements of uniformity and 

independence.

0.0000

0.0878

0.1756

0.2635

0.3513

0.4391

0.000

0.010

0.020

0.030

0.040

0.0

55.2

110.4

165.6

220.8

276.0

331.2

386.4

441.6

496.8

552.0

607.3

662.5

717.7

Peak p (normalized) 

Relative Frequency 

Duration d (ms) 

Empirical Probability Distribution Function 

Figure 12: Empirical Probability Density Function

12
0.00

0.02

0.04

0.06

0.08

0.10

0.12

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 More

Relative Frequency

Range

PRNG Probability Distribution Function

0.00

0.10

0.20

0.30

0.40

0.50

0.60

0.70

0.80

0.90

1.00

Relative Frequency

Range

PRNG Cumulative Distribution Function

Figure 14: Probability distribution function for sample output of pseudo-random number generator

Figure 13: Cumulative distribution function for sample of pseudo-random number generator

13
4.3. Overstroke Sample Generation

The experimental samples and the samples generated by the R ejection Sampling and Inverse 

Transform Sampling techniques are shown in Figures 16 to 18. The corresponding marginal CDFs for the 

overstroke peak p and duration d are given in Figures 19 and 20, respectively. Visual inspection of these 

figures appears to indicate that the Inverse Transform Sampling technique has slightly better fidelity. The 

following section gives objective measures of goodness-of-fit for both simulation techniques.

0.00

0.10

0.20

0.30

0.40

0.50

0.60

0.70

0.80

0.90

1.00

0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00

PRNG2

PRNG1

PRNG Scatter Plot

Figure 15: Scatter plot for two PRNGs

14
0.000

0.100

0.200

0.300

0.400

0.500

0.600

0.0 100.0 200.0 300.0 400.0 500.0 600.0 700.0 800.0

p

d (ms)

Experimental Sample

0.000

0.100

0.200

0.300

0.400

0.500

0.600

0.0 100.0 200.0 300.0 400.0 500.0 600.0 700.0 800.0

p

d (ms)

Rejection Sampling

Figure 16: Experimental overstroke samples

Figure 17: Samples generated by Rejection Sampling technique

15
0.000

0.100

0.200

0.300

0.400

0.500

0.600

0.0 100.0 200.0 300.0 400.0 500.0 600.0 700.0 800.0

p

d (ms)

Inverse Transform Sampling

0.00

0.20

0.40

0.60

0.80

1.00

1.20

0.0000.0370.0730.1100.1460.1830.2200.2560.2930.3290.3660.4030.4390.476

Cumulative Relative Frequency

p

Cumulative Distributions of Overstroke Peak p

Experimental

Rejection

Sampling

Inverse PIT

Figure 19: Marginal CDF of overstroke peak p for experimental and simulation samples

Figure 18: Samples generated by Inverse Transform Sampling technique

16
4.4. Goodness-of-Fit Measures

The goodness of fit (GoF) of the simulation samples was evaluated quantitatively using the ChiSquare test [16, 17] applied to the marginal CDFs for p and d of the generated samples using 10 cells, 

corresponding to 9 degrees of freedom. The expected observed frequency was based on the experimental 

sample. The cell (or bin) widths were defined such that each bin covered 10% of the points in the original 

experimental sample. A s the Chi-Square statistics in Table 1 show, the Inverse Transform Sampling 

technique performed much better for both p and d than the Rejection Sampling technique. It is not clear 

why Inverse Transform Sampling had better results, though we suspect this may be due t o some sort of 

smoothing effect in the Inverse Probability Integral Transform (IPIT) that minimizes the slight residual 

clustering patterns (i.e., correlation) in the PRNG outputs visible in Figure 15. The Rejection Sampling 

technique simply accepts or re jects generated (p, d) poin ts with no significant mathematical 

transformation between the PRNG outputs and the points.

Table 1: Chi-Square Test results for marginal CDFs (9 degrees of freedom)

Overstroke Parameter Simulation Technique Chi-Square Statistic Significance

Peak p Rejection Sampling 47.81 2.78e-07

Inverse Transform Sampling 7.36 0.60

Duration d Rejection Sampling 376.22 1.51e-75

Inverse Transform Sampling 15.56 0.08

0.00

0.20

0.40

0.60

0.80

1.00

1.20

0.0

55.2

110.4165.6220.8276.0331.2386.4441.6496.8552.0607.3662.5717.7

Cumulative Relative Frequency

d

Cumulative Distributions of Overstroke Duration d

Experimental

Rejection

Sampling

Inverse PIT

Figure 20: Marginal CDF of overstroke duration d for experimental and simulation samples

17
5. Conclusions

We have described a general modeling and simulation approach for a typical disturbance pattern that a

digital device may experience in an electromagnetic reverberation chamber. The modeling approach is 

based on a three-dimensional histogram using bilinear interpolation to build a continuous surface for the 

empirical probability distribution function. We cons idered two simulation approaches extended to the 

case of a bivariate distribution: Rejection Sampling and Inverse Transform Sampling. The results showed 

that the I nverse Transform Sampling technique can generate samples that most closely match the 

empirical distribution. 

We intend to apply this modeling and simulation approach to generate random disturbance patterns as 

stimulus in a sy stem fault-response simulator. The purpose of the simulator is to gain a th orough 

understanding of the relation between internal upsets and the propagated effects observable at the external 

system interface. To achieve this, we needed to develop a technique to map disturbances to internal 

perturbations of the sy stem failure units (see Figure 3). We a lso need to d evelop error propagation 

models that de termine the system-structure-dependent transformation from internal perturbations to 

external disruptions (see Figure 2). The ultimate goal is to develop a model-based methodology to assess 

system fault response and guide design efforts for enhanced resilience in complex safety-critical 

distributed systems.

References

[1] Federal Aviation Administration. http://www.faa.gov/nextgen/ Accessed February 24, 2014.

[2] Driscoll, Kevin, et al.: Data Network Criteria Evaluation Report. Fed eral Aviation Administration, 

DOT/FAA/AR-09/27, July 2009.

[3] Torres-Pomales, Wilfredo: An Approach for the Assessment of System Upset Resilience. NASA/TM2013-217798, January 2013.

[4] Belcastro, Celeste M.: Laboratory Test Methodology for Evaluating the Effects of Electromagnetic 

Disturbances on Fault-Tolerant Control Systems. NASA Technical Memorandum 101665, 

November 1989.

[5] Fuller, Gerald L.: Understanding HIRF – High Intensity Radiated Fields. Aviation Communications, 

Inc., Leesburg, VA, 1995, p. 7-2.

[6] Hess, Richard: Computing Platform Architectures for Robust Operation in the Presence of Lightning 

and Other Electromagnetic Threats. 16th Digital Avionics Systems Conference (DASC), 1997.

[7] Torres-Pomales, W.: Analysis of the Radiated Field in an Electromagnetic Reverberation Chamber as 

an Upset-Inducing Stimulus for Digital Systems. NASA/TM-2012-217789, December 2012.

[8] Williams, R. A.: The NASA High Intensity Radiated Fields Laboratory. 16th Digital Avionics Systems 

Conference, October 1997.

[9] Torres-Pomales, Wilfredo: Characterization of the HIRF Susceptibility Threshold for a Prototype 

Implementation of an Onboard Data Network. NASA/TM-2012-217754, August 2012.

18
[10] Yates, A.M.; et al.: High-Intensity Radiated Field Fault-Injection Experiment for a Fault-Tolerant 

Distributed Communication System. IEEE/AIAA 29th Digital Avionics Systems Conference (DASC), 

October 2010, pp. 4.E.301 – 4.E.3-15.

[11] Mendenhall, W. and Sincich, T.: Statistics for Engineering and the Sciences. 5th Edition, P rentice 

Hall, 2006.

[12] Nelsen, Roger B.: An Introduction to Copulas. 2nd Edition, Springer, 2007.

[13] van der Vaart, A. W.: Asymptotic Statistics. Cambridge University Press, 2000.

[14] Scott, David W.: Multivariate Density Estimation: Theory, Practice, and Visualization. 1 st

Edition, Wiley, 1992.

[15] Fishman, George: Monte Carlo: Concepts, Algorithms, and Applications. Springer, 2011.

[16] Barnes, J. Wesley: Statistical Analysis For Engineers. Prentice Hall, 1988.

[17] Greenwood, Priscilla E.; and Nikulin, Michael S. A Guide to Chi-Squared Testing. Wiley-Interscience, 

1996.

19
Abbreviations

CDF Cumulative Distribution Function

eCDF Empirical Cumulative Distribution Function

ePDF Empirical Probability Distribution Function

GoF Goodness of Fit

HIRF High Intensity Radiated Field

IPIT Inverse Probability Integral Transform

ms Millisecond

NASA National Aeronautics and Space Administration

PD Probability Distribution

PDF Probability Distribution Function

PIT Probability Integral Transform

PRNG Pseudo-Random Number Generator

V&V Validation and Verification

20
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Memorandum

 4. TITLE AND SUBTITLE

Modeling and Simulation of Upset-Inducing Disturbances for Digital 

Systems in an Electromagnetic Reverberation Chamber 

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Torres-Pomales, Wilfredo

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

L-20380

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 33

Availability: NASA CASI (443) 757-5802

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

This report describes a modeling and simulation approach for disturbance patterns representative of the environment 

experienced by a digital system in an electromagnetic reverberation chamber. The disturbance is modeled by a multi-variate 

statistical distribution based on empirical observations. Extended versions of the Rejection Samping and Inverse Transform 

Sampling techniques are developed to generate multi-variate random samples of the disturbance. The results show that Inverse 

Transform Sampling returns samples with higher fidelity relative to the empirical distribution. This work is part of an ongoing

effort to develop a resilience assessment methodology for complex safety-critical distributed systems. 

15. SUBJECT TERMS

Modeling; Safety; Simulation

18. NUMBER

 OF 

 PAGES

27

19b. TELEPHONE NUMBER (Include area code)

(443) 757-5802

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

 534723.02.02.07.30 

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/TM-2014-218243

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 04 - 2014