Structured Random Field Extender (SRFE)
From BoimWiki
This is an algorithm for generating a 2D random field, as displayed to the right, which adheres to a given structure function. The primary intended purpose for this procedure is to generate Kolmogorov phase screens, which represent index of refraction variations in an atmosphere described by Kolmogorov turbulence.
Contents |
History
It is based primarily off of the paper "Method for simulating infinitely long and non stationary phase screens with optimized memory storage" by François Assémat and Richard W. Wilson with Eric Gendron which appeared in Optics Express Vol. 14, No. 3, 6 February 2006, page 988. I have also had some conversations with Dr. Tim Clark, who later published a paper, perhaps in the same journal with Dr. David Fried with an extension on the topic. AFAIK, their contribution was to extend the prediction input point back farther from the front row, and introduce the idea of sampling a subset of the field for prediction, since using all available data would take too much time and memory.
In both of the papers above, they predict an entire new row of the random field at a time. The primary extension of the algoritm presented here is that it extends the field a single pixel at a time. This simplifies the math significantly. Many vector and matrix operations become scalar and vector operations respectively. Once a single pixel can be predicted, the same procedure can be repeated for every pixel in a row to be expanded. This allows for smaller prediction sample fields (stencils), lower memory requirements, and faster operation.
Derivation
Assume we have an answer of the form Failed to parse (unknown error): Z^{T} A + B \beta = x . Where Z are the values of previously known samples of the random field, at known locations. A is a prediction vector for the mean of the new sample, x, and B is the standard deviation of the new sample X, and Failed to parse (unknown error): \beta
is an RV of mean 0, standard deviation 1.
Pre-multiply by Z: Failed to parse (unknown error): Z Z^T A + Z B \beta = Zx
Take the expected value,
Failed to parse (unknown error): <Z Z^T>A + B<Z\beta> = <Zx>
Note that Failed to parse (unknown error): \beta
and Z are uncorrelated, expected value 0. Failed to parse (unknown error): <Z Z^T>A = <Zx>
The correlation matrix for Z, and the correlation of Z to the new pixek, x are known (we are given a structure function).
We can solve for A.
If we return to the original solution assumption and square it, Failed to parse (unknown error): x^2 = A^T Z Z^T A + 2 A^T Z \beta B + B^2 \beta^2
Take the expected value. Since Failed to parse (unknown error): \beta
and Z are uncorrelated, the middle term drops out, the expected value of Failed to parse (unknown error): \beta^2
is one, hence Failed to parse (unknown error): <xx> = A^T <Z Z^T> A + B^2
Plug in the previous solution for the second A, the ZZ ZZ-1 cancel, and we get:
Failed to parse (unknown error): B^2 = <xx> - A^T <xZ>
Use
The primary use for these Kolmogorov phase screens is in wave-optics propagation simulations. The ability to expand a large phase screen indefinitely allows for consistent simulations which involve a slew across large swaths of atmosphere.
The technique can be generalized for filling in sparse random fields, or expanding a random field in any direction. The row-wise extension is provided as a simple example. Code in the lower-level interface could be re-used for arbitrary expansion procedures.
The structure function used for extension is itself a parameter. The only structure function provided is Kolomogorov, but others could be readily implemented. The user must provide a covariance function, which is a function of a 2D offset vector. In the case of Kolmogorov, the covariance is a function of distance, and has no directional dependence. However, covariance functions which are direction dependent could be implemented.
Implementation
Source Code in C is provided to demonstrate the procedure. C was chosen as a fairly least-common-denominator language which can be called by many higher level languages. It has a dependency on LAPACK to perform a matrix inversion.
Row Extension Interface
See the Command Line section below for a description of how to generate coefficients for a row extender.
The program in writePhaseScreen.c provides an example of how to extend rows of a Kolmogorov phase screen. It consists of three simple calls:
SRFEreLoad(nCol, dRow, coefFileNameFormat, &re); SRFEreExtend(&re); /* extend field by 1 row. */ SRFEreDelete(&re); /* free up resources for the row extender */
The SRFEreLoad() call loads the row extension coefficients from files generated by the KolRowGen program. The SRFEreExtend() call adds one new row to the ring buffer stored in the random field extender. Call SRGEreDelete() to free the resources used for the extender.
It can take several minutes to derive the predictor coefficients (KolRowGen). However, when running the extender, the time necessary to generate a new row is negligible.
Following is a code snippet, similar to what is seen in writePhaseScreen.c, which uses this interface:
#include "SRFE.H"
int dRow; /* number of rows to look back for extending phase screen */
int nCol; /* width of phase screen to generate */
const char *fNamFmt;
SRFErowExtender re;
...
double kolParam = dRow;
SRFEreLoad(nCol, dRow, fNamFmt, &re);
...
SRFEreEextend(&re); /* extend field by 1 row. oldest row is lost from buffer */
...
/* re.row[0] is a float*, newest row on phase screen.
re.row[1] is the previous row ...
re.row[re.dRow] is the oldest row in the buffer */
...
SRFEreDelete(&re); /* free up resources for the row extender */
--place descriptions of lower level calls in the API here--
Contact 
for more information, or help with using this code.
If you wish to publish a paper based on this work, contact me (Aaron).
I will be happy to share authorship of a paper with anybody who is willing
to do some additional validation on the algorithm, like verification
that the random field has the desired structure function.
Prediction Stencil
| This animation shows the points used for prediction of each pixel on the new row.
The location of each pixel used for the Z vector is plotted, with the color of each plotted pixel corresponding to the associated A coefficient magnitude. Note that most of the coefficients for these predictors are near zero, except for close to the point. This makes me wonder of such a large or dense prediction stencil is really necessary. I believe that Clark has shown that the large extent of pixels examined for prediction helps to enforce the structure function over larger distances. For some reason, transitivity of correlations is not enough to uphold the structure function over large distances. His measurements indicate that the structure function rolls off where points are not used in the predictor. I am not sure this would be necessary for many applications. Experiments should be done with more sparse prediction arrays. It is hard for me to believe that all those near-zero coefficients have much effect. |
 |
 |
Command Line Interface
The prediction coefficient generator has been de-coupled from phase-screen extension. This is a good thing since coefficient generation is a lengthy process, and most users will want to do this once, then re-use the predictors thousands of times. It also removes some significant dependencies from the phase-screen generation routines.
It is recommended that you create a folder for the hundreds (thousands) of small coefficient files. The width and depth of the Region of Support are also needed for future prediction runs, so it is suggested that these numbers be encoded into the folder name. Here are the commands for an example coefficient generation run:
mkdir SRFE1024_2047 KolRowGen 1024 2047 100 SRFE1024_2047/pred%04d.dat
Although the code to write a long phase screen was intended as an example of how to generate rows of a phase screen, in a circular buffer, for extended use, it may be of utility to some users. This would be the command to generate a 1024x8196 pixel phase screen
writePhaseScreen 1024 SRFE1024_2047/pred*04d.dat ps1024x8196.bin 8196 4000
Note that the same prediction coefficient file format argument was used here as was used for generation. Note also that the same width (1024) parameter was also used. The last argument is the number of lines to generate, but not save. It is a good idea to cycle through at least one or two buffers worth of data (2047 here) to ensure that the statistics have settled, and estimates are not biased by the initial condition of an all zero field. 4000 was chosen here since it is roughly two times the buffer depth of 2047.
Future Work
Some parts of the work are not fully generalized. There was some normalization of the Kolmogorov phase screen which might be unnecessary. I need to do some more testing to see if this dependency can be eliminated. I am also eager to hear about other structure functions that might be of interest.
The random part of these fields are generated with a Gaussian RV. In some discussions with Clark, I have reason to believe that a double-exponential function might be a more appropriate distribution. This could easily be put into the algoruthm. The only requirement is that the random part have the indicated Standard Deviation (second moment). The function which has this second moment could be double exponential as easily as Normal.
To farther demonstrate and validate the procedure, code for another common random field structure, Gaussian was implemented. The procedure seemed to work for some Gaussian structure function parameters, but not others. Some Gaussian based correlation matrices were not invertable. Others were singular, but seemed to generate non-stable predictors. I have experimented with changes to the prediction stencil, thinking that the system may have been over-determined for the faster roll-off Gaussian structure function. These experiments are on-going, and not yet totally successful. Completing the test of a Gaussian random field should be a good validation of the procedure (and it's limitations)
Version History
- 091103 initial draft. integrated extension and generation
- 091104 separated generation code from extension code
