# Exploring the Effectiveness of Self-supervised Learning and Classifier
Chains in Emotion Recognition of Nonverbal Vocalizations

Detai Xin  Shinnosuke Takamichi  Hiroshi Saruwatari

###### Abstract

We present an emotion recognition system for nonverbal vocalizations (NVs)
submitted to the ExVo Few-Shot track of the ICML Expressive Vocalizations
Competition 2022. The proposed method uses self-supervised learning (SSL)
models to extract features from NVs and uses a classifier chain to model the
label dependency between emotions. Experimental results demonstrate that the
proposed method can significantly improve the performance of this task
compared to several baseline methods. Our proposed method obtained a mean
concordance correlation coefficient (CCC) of $0.725$ in the validation set and
$0.739$ in the test set, while the best baseline method only obtained $0.554$
in the validation set. We publicate our code at https://github.com/Aria-K-
Alethia/ExVo to help others to reproduce our experimental results.

Nonverbal Vocalizations, Emotion Recognition

##  1 Introduction

Nonverbal vocalizations (NVs), also called affect bursts, refer to short and
expressive vocalizations containing no linguistic information like laughter,
sobs, and screams (Scherer, 1994; Trouvain & Truong, 2012). NVs are important
for spoken language processing, since (1) NVs play an important role in
expressing emotions in spoken languages (Hall et al., 2009; Scherer & Scherer,
2011), and (2) they are common components of human communication existing in
different cultures and languages (Sauter et al., 2010). Although both
emotional prosody and NVs contribute to emotional expressions in speech, NVs
are ignored by most previous research on speech emotions (Lima et al., 2013),
which necessitates further work in this field.

The ICML Expressive Vocalizations Competition (hereafter ExVo) launched in
2022 aims to develop technologies for the recognition, generation, and
personalization of NVs (Baird et al., 2022). The competition includes three
tracks: ExVo Multi-Task, ExVo Generate, and ExVo Few-Shot, which correspond to
the three goals of recognizing, generating, and personalizing NVs,
respectively. Specifically, the ExVo Multi-Task track aims to recognize not
only emotions but also demographic information like the age and native country
of the speakers from NVs. In the ExVo Generate track, the participants need to
generate NVs of various emotions. Similar to the Multi-Task track, the ExVo
Few-Shot track also aims to recognize emotions from NVs, but further requires
the prediction systems to adapt to new speakers which are not in the training
set. All tracks are evaluated by appropriate subjective and objective metrics
to reflect the performance of the proposed systems.

In this paper, we describe our emotion recognition system for NVs submitted to
the ExVo Few-Shot track. Our system consists of two components: a feature
extractor and a classifier chain (CC). As the feature extractor, we use a
self-supervised learning (SSL) model like Wav2vec2 (Schneider et al., 2019;
Baevski et al., 2020) or HuBERT (Hsu et al., 2021a). The extracted feature is
then fed to a classifier chain, which predicts the score for each emotion
sequentially. The proposed classifier chain predicts the score by conditioning
on both the extracted feature and the predicted scores in previous steps, so
can utilize information from label dependency. We conducted comprehensive
experiments to verify the effectiveness of these two components, and show that
the proposed method can significantly improve the performance compared to
several baseline methods. Our contributions can be summarized as follows:

  * •

We propose an emotion recognition system for NVs using SSL models and CC, and
conduct experiments to show the effectiveness of the proposed method.

  * •

We conduct experiments to show the best SSL models for this task.

  * •

We analyze the prediction results and give insights for future research.

##  2 The ExVo Few-Shot track

The task of the ExVo Few-Shot track is recognizing emotions from NVs for
unseen speakers that are not in the training set. During the training phase,
the participants can train their emotion recognition models on the data of the
seen speakers. In the test phase, two samples for each unseen speaker will be
provided by the organizer, and the participants can then adapt their models to
unseen speakers by few-shot learning.

ExVo provides a large-scale multilingual (Chinese, English, Spanish) NVs
dataset (Cowen et al., 2022), which contains approximately $36$ hours NVs
uttered by $1,702$ speakers from the USA, China, South Africa, and Venezuela.
Ten emotions consisting of amusement, awe, awkwardness, distress, excitement,
fear, horror, sadness, surprise, and triumph are covered by the dataset. Each
NV is annotated with ten emotion scores ranging from $0$ to $1$ by crowd-
sourcing. Some statistics of the dataset are summarized in Table 1.

Table 1: Statistics of the dataset. Some terms are not available for the participants. | Train | Validation | Test | $\sum$
---|---|---|---|---
Duration (hrs) | $12.32$ | $12.10$ | $12.37$ | $36.78$
Speakers | $571$ | $568$ | $563$ | $1702$
Audio clips | $19,990$ | $19,396$ | $19,815$ | $59201$
Audio clips per speaker |  $4$-$201$ (Avg. $35$) |  $4$-$168$ (Avg. $34$) | – | –

The participants should submit the predictions of their models on test
samples. Concordance correlation coefficient (CCC) is used to evaluate the
performance, which can measure the agreement between the ground-truth (GT) and
predicted emotion scores.

##  3 Proposed method

In this section, we describe the proposed method. We first introduce the
general architecture of the proposed method and then describe each component
separately.

###  3.1 Architecture

The architecture of the proposed method is illustrated in Figure 1.
Specifically, the NVs are first fed to the SSL model to extract sequential
features from them. Then, an attentive pooling module is used to collect
information over the time axis and convert the sequential features into fixed-
length features. Finally, a CC is used to predict the score for each emotion
sequentially. Each emotion has a separate predictor, which contains a fully
connected (FC) layer followed by a sigmoid activation. The predictor in the CC
not only uses the extracted features but also the predicted scores of previous
classifiers as input, which can thus utilize information of label dependency.

Figure 1: Architecture of the proposed method. FC denotes a fully connected
layer.

###  3.2 Self-supervised learning models

SSL is a popular and powerful method to leverage large-scale unlabeled data.
The idea of SSL is to set a sophisticated pretext task to learn nontrivial
data representations. Recently, several works based on SSL have been proposed
in speech processing, which showed promising results in speech-to-text and
other speech-based tasks (Schneider et al., 2019; Baevski et al., 2020;
Conneau et al., 2020; Hsu et al., 2021b, a). Although previous work has shown
that pretrained SSL models were effective on speech emotion recognition in
verbal communication(Pepino et al., 2021; Chen & Rudnicky, 2021; Wang et al.,
2021; Wagner et al., 2022), since most of the SSL models were trained on
speech corpora that rarely contain NVs, whether these models are effective on
NVs remains unknown. Therefore, in the experiments, we select several SSL
models and verify their effectiveness in the task of recognizing emotions from
NVs.

###  3.3 Classifier chain

The task of the ExVo Few-Shot track is a multivariate regression problem.
While similar tasks like continuous speech emotion recognition usually use
separate predictors to predict the emotions independently (Atmaja & Akagi,
2021), it is noteworthy that the emotion labels in the provided dataset have
an intrinsic dependency. For example, if the score of amusement is high, the
scores of negative emotions like sadness and fear may be low, which implies
the possibility to utilize information from label dependency. To this end, we
propose to use CC (Read et al., 2011), which has been proved to be useful for
modeling label dependency in multi-label classification (Dembszynski et al.,
2010; Dembczynski et al., 2010; Nam et al., 2017).

Formally, denoting the predictor of the $i$-th emotion and the extracted
feature as $f_{i}(\cdot)$ and $z$, respectively, the emotion score
$\widehat{y}_{i}$ is computed by:
$\widehat{y}_{i}=\sigma(f_{i}(z\oplus\widehat{y}_{<i})),$ where
$\sigma(\cdot)$ is the sigmoid function, $\oplus$ is the vector concatenation
operator, and $\widehat{y}_{<i}$ indicates a vector concatenating previously
predicted emotion scores. During training, we feed GT emotion scores $y_{<i}$
to the predictor to prevent error propagation.

While it is possible to use powerful sequential models like a recurrent neural
network to implement the predictor $f_{i}$ in CC (Nam et al., 2017), in the
preliminary experiments we found simple linear CC was more stable than other
models, thus we use linear CC in the proposed method. Chain order is another
critical problem for CC since it determines what information can be utilized
by each classifier (Read et al., 2021). In our implementation, we adopt a
heuristic method in which we first train ten separate base predictors for each
emotion, then use the descending order of the performance of the base
predictors as the chain order. This method can intuitively alleviate the error
propagation problem by first predicting emotions with high confidence, hence
can improve the performance (Read et al., 2021).

###  3.4 Loss function

We use CCC as the objective for the model training. Previous work has
demonstrated that CCC was better for speech emotion recognition than error-
based loss functions like L1 loss (Atmaja & Akagi, 2021). Generally, given a
set of paired data ($\\{s^{j}\\}_{j=1}^{n},\\{t^{j}\\}_{j=1}^{n}$) with length
$n$, the CCC between them is defined as:
$\mathrm{CCC}(\\{s^{j}\\}_{j=1}^{n},\\{t^{j}\\}_{j=1}^{n})=\frac{2\sigma^{2}_{st}}{\sigma^{2}_{s}+\sigma^{2}_{t}+(\mu_{s}-\mu_{t})^{2}}.$
In the proposed method the objective value of a mini-batch is computed by
averaging CCC values of all emotions. Formally, denoting the number of
emotions and the mini-batch size as $C$ and $B$, respectively, the loss
function of the proposed model is defined as:
$\mathcal{L}=\frac{1}{C}\sum_{i=1}^{C}\mathrm{CCC}(\\{y_{i}^{j}\\}_{j=1}^{B},\\{\widehat{y}_{i}^{j}\\}_{j=1}^{B}),$
where $y_{i}^{j}$ and $\widehat{y}_{i}^{j}$ are GT and predicted scores of the
$i$-th emotion of the $j$-th sample in the mini-batch, respectively.

###  3.5 Data augmentation

To improve the robustness of the model, we additionally use data augmentation.
Specifically, we use two strategies: pitch-shifting and speaking-rate-changing
(Saeki et al., 2022). Pitch-shifting raises or lowers the pitch of NVs, and
can change the speaker identity of NVs. Speaking-rate-changing slows down or
speeds up the NVs. We tune the shifting range and the speaking-rate-changing
range so that the emotion of the augmented NVs has little difference from the
original ones.

##  4 Experiments

Table 2: Average CCC of different SSL models in $5$-fold cross-validation. Bold indicates the best score. SSL Model | CCC
---|---
Wav2vec2-base | $0.699\pm 0.013$
HuBERT-base | $0.703\pm 0.013$
Wav2vec2-robust | $0.709\pm 0.012$
Wav2vec2-large | $0.709\pm 0.009$
HuBERT-large | $0.718\pm 0.011$
XLSR | $\mathbf{0.722\pm 0.012}$
Table 3: Average CCC of the proposed method using CC and data augmentation in 5-fold cross-validation. Bold indicates the best score. Model | CCC
---|---
CC HuBERT-large | $0.722\pm 0.012$
CC HuBERT-large Aug. | $0.722\pm 0.012$
CC XLSR | $0.724\pm 0.012$
CC XLSR Aug. | $\mathbf{0.726\pm 0.013}$
Table 4: CCC of the proposed and baseline methods for each emotion on all validation data. Bold indicates the best score with p-value $<0.05$. Model | Awe | Excite. | Amuse. | Awkward. | Fear | Horror | Distress | Triumph | Sadness | Surprise | Average
---|---|---|---|---|---|---|---|---|---|---|---
DeepSpectrum | $0.575$ | $0.424$ | $0.512$ | $0.318$ | $0.497$ | $0.480$ | $0.356$ | $0.263$ | $0.412$ | $0.621$ | $0.446$
eGeMAPS | $0.592$ | $0.420$ | $0.525$ | $0.302$ | $0.527$ | $0.511$ | $0.365$ | $0.263$ | $0.407$ | $0.604$ | $0.452$
BoAW | $0.639$ | $0.470$ | $0.583$ | $0.383$ | $0.566$ | $0.548$ | $0.457$ | $0.347$ | $0.472$ | $0.659$ | $0.512$
ComParE | $0.680$ | $0.503$ | $0.638$ | $0.403$ | $0.604$ | $0.587$ | $0.494$ | $0.411$ | $0.526$ | $0.704$ | $0.554$
CC XLSR Aug. | $\mathbf{0.801}$ | $\mathbf{0.687}$ | $\mathbf{0.796}$ | $\mathbf{0.608}$ | $\mathbf{0.755}$ | $\mathbf{0.735}$ | $\mathbf{0.685}$ | $\mathbf{0.674}$ | $\mathbf{0.703}$ | $\mathbf{0.804}$ | $\mathbf{0.725}$

###  4.1 Experimental Setup

We constructed several baseline systems based on the features provided by the
organizer (Baird et al., 2022). These features include: the $6373$-dimensional
ComParE set from the 2016 Computational Paralinguistics Challenge (Schuller et
al., 2016) and the $88$-dimensional extended Geneva Minimalistic Acoustic
Parameter Set (eGeMAPS) (Eyben et al., 2015) extracted by openSMILE toolkit
(Eyben et al., 2010), the Bag-of-Audio-Words (BoAW) representations with a
$2000$ codebook size extracted by openXBOW toolkit (Schmitt & Schuller, 2017)
from the low-level descriptors of ComParE set, the $4096$-dimensional spectrum
representations extracted by the DeepSpectrum toolkit (Amiriparian et al.,
2017). We used multi-layer neural networks with LeakyReLU activation (Maas et
al., 2013) to process each of the features. Batch normalization (Ioffe &
Szegedy, 2015) was used in each layer to normalize the features.

To verify the effectiveness of SSL models and find the best model for NVs, we
selected six SSL models: Wav2vec2-{base, large} (Baevski et al., 2020),
HuBERT-{base, large} (Hsu et al., 2021a), Wav2vec2-robust (Hsu et al., 2021b),
and XLSR (Conneau et al., 2020). Wav2vec2-{base, large} and HuBERT-{base,
large} are common models used in previous work (Chen & Rudnicky, 2021; Wang et
al., 2021). The Wav2vec2-{base, large} and HuBERT-large were trained on the
Librispeech (Panayotov et al., 2015) corpus containing $960$ hours of audio,
while the HuBERT-large model was trained on the Libri-Light (Kahn et al.,
2020) corpus containing $60$k hours of audio. These corpora are all from clean
audiobooks that rarely contain NVs. Also, since the NVs in the dataset were
recorded by the speakers themselves and sometimes contain noises, we selected
Wav2vec2-robust that was trained on not only the Libri-Light corpus but also
noisy telephone corpora including Switchboard (Godfrey & Holliman, 1997),
Fisher (Cieri et al., 2004), and CommonVoice (Ardila et al., 2019). Finally,
to handle multilingual NVs we select XLSR, which is a multilingual version of
Wav2vec2 trained on $56$k hours of multilingual audio (multilingual
Librispeech (Pratap et al., 2020), CommonVoice, and Babel (Gales et al.,
2014)) covering $53$ languages.

Instead of using the original train-validation split, we first combined the
train and validation sets and then used $5$-fold cross-validation to train all
models. To ensure the validation results can reflect the generalization
ability of the models, we split the data by speaker so that the speakers of
the validation set of each split are unseen for the training set. Each split
had about $912$ seen speakers and $227$ unseen speakers. All silence in the
audios was trimmed by voice activity detection. The chain order was set to
awe, surprise, amusement, fear, horror, sadness, distress, excitement,
triumph, and awkwardness. We used Adam (Kingma & Ba, 2015) to optimize the
model, the learning rate was set to $1\mathrm{e}{-4}$ for the CC and
$1\mathrm{e}{-5}$ for the feature extractor. The learning rate was halved once
the loss value on the validation set did not decrease. The batch size was set
to $16$ when Wav2vec-base or HuBERT-base was used, otherwise, it was set to
$8$. We used WavAugment (Kharitonov et al., 2020) to implement the data
augmentation, the pitch-shifting range was set to $[-300,300]$ cents, and the
speaking-rate-changing range was set to $[0.8,1.2]$. The minimal and maximal
number of epochs were set to $10$ and $50$, respectively. Early stopping was
used to prevent overfitting with a patience of $10$ epoch. We selected the
model that had the best CCC on the validation set. To adapt the models to
unseen speakers, we further fine-tuned the models on the $2$ samples of each
test speaker, which roughly took $10$ epochs for convergence. The learning
rate was set to $1\mathrm{e}{-6}$ for the whole network.

###  4.2 Evaluations of SSL models

We first evaluate all SSL models to verify their effectiveness and find the
best model for NVs. To avoid the the influence on performance from the CC, in
this experiment the CC is replaced by a single FC layer. The results are shown
in Table 2. It can be seen that the XLSR model obtained the best performance.
We suppose this is because the multilingual XLSR model can handle more
phonetic tokens of different languages, unlike other SSL models that are only
trained on English corpora. Also, large models are always better than base
models, which is consistent with the results of previous work (Wagner et al.,
2022).

###  4.3 Evaluations of CC and data augmentation

We select the top-2 SSL models (XLSR and HuBERT-large) in the previous
experiment as the feature extractors to evaluate CC and data augmentation
(“Aug.”). The results are shown in Table 3. It can be seen that CC
consistently improves the performance compared to the results in Table 2 The
data augmentation brings improvements for the XLSR model but failed to improve
the performance of the CC HuBERT-large model.

###  4.4 Evaluations for each emotion

We then evaluate the proposed method for each emotion. We selected the best
performing model CC XLSR Aug. obtained in previous experiments. All baseline
models were trained using cross-validation for comparison. After training, we
combined the inference results on all $5$ splits of each model and computed
the CCC value for each emotion. The result is demonstrated in Table 4. It can
be seen that the proposed method significantly outperformed all baseline
methods on all emotions, which indicates the effectiveness of SSL models and
CC. Besides, it can be observed that all models have difficulties to recognize
some emotions like awkwardness, triumph, distress.

###  4.5 Evaluations on test set

We finally fine-tuned the best model (CC XLSR Aug.) on the test samples
provided by the organizer. Our best CCC on the test set is $0.739$, which
demonstrates the proposed method has a strong generalization ability.

##  5 Conclusions

This paper described an emotion recognition system for NVs submitted to the
Few-Shot track of the ICML ExVo competition. The proposed method uses a SSL
model to extract contextual speech representations from NVs, and uses a CC to
predict emotion scores by utilizing the features and the emotion scores
predicted in previous steps together. Experimental results demonstrated that
the proposed method significantly outperformed several baseline methods.

Acknowledgements: This work was supported by JST SPRING, Grant Number
JPMJSP2108. Part of this work was also supported by JSPS KAKENHI Grant Number
21H04900 (for implementation) and JST Moonshot R&D Grant Number JPMJPS2011
(for evaluation).

## References

  * Amiriparian et al. (2017) Amiriparian, S., Gerczuk, M., Ottl, S., Cummins, N., Freitag, M., Pugachevskiy, S., Baird, A., and Schuller, B.  Snore sound classification using image-based deep spectrum features.  _Proc. Interspeech 2017_ , pp. 3512–3516, 2017.
  * Ardila et al. (2019) Ardila, R., Branson, M., Davis, K., Henretty, M., Kohler, M., Meyer, J., Morais, R., Saunders, L., Tyers, F. M., and Weber, G.  Common voice: A massively-multilingual speech corpus.  _arXiv preprint arXiv:1912.06670_ , 2019.
  * Atmaja & Akagi (2021) Atmaja, B. T. and Akagi, M.  Evaluation of error-and correlation-based loss functions for multitask learning dimensional speech emotion recognition.  In _Journal of Physics: Conference Series_ , volume 1896, pp. 012004. IOP Publishing, 2021.
  * Baevski et al. (2020) Baevski, A., Zhou, Y., Mohamed, A., and Auli, M.  wav2vec 2.0: A framework for self-supervised learning of speech representations.  _Advances in Neural Information Processing Systems_ , 33:12449–12460, 2020.
  * Baird et al. (2022) Baird, A., Tzirakis, P., Gidel, G., Jiralerspong, M., Muller, E. B., Mathewson, K., Schuller, B., Cambria, E., Keltner, D., and Cowen, A.  The icml 2022 expressive vocalizations workshop and competition: Recognizing, generating, and personalizing vocal bursts.  _arXiv preprint arXiv:2205.01780_ , 2022.
  * Chen & Rudnicky (2021) Chen, L.-W. and Rudnicky, A.  Exploring wav2vec 2.0 fine-tuning for improved speech emotion recognition.  _arXiv preprint arXiv:2110.06309_ , 2021.
  * Cieri et al. (2004) Cieri, C., Miller, D., and Walker, K.  The fisher corpus: A resource for the next generations of speech-to-text.  In _LREC_ , volume 4, pp. 69–71, 2004.
  * Conneau et al. (2020) Conneau, A., Khandelwal, K., Goyal, N., Chaudhary, V., Wenzek, G., Guzmán, F., Grave, É., Ott, M., Zettlemoyer, L., and Stoyanov, V.  Unsupervised cross-lingual representation learning at scale.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pp. 8440–8451, 2020.
  * Cowen et al. (2022) Cowen, A., Baird, A., Tzirakis, P., Opara, M., Kim, L., Brooks, J., and Metrick, J.  The hume vocal burst competition dataset (h-vb) — raw data [exvo: updated 02.28.22] [data set].  _Zenodo_ , 2022.  doi: https://doi.org/10.5281/zenodo.6308780.
  * Dembczynski et al. (2010) Dembczynski, K., Cheng, W., and Hüllermeier, E.  Bayes optimal multilabel classification via probabilistic classifier chains.  In _ICML_ , 2010.
  * Dembszynski et al. (2010) Dembszynski, K., Waegeman, W., Cheng, W., and Hüllermeier, E.  On label dependence in multilabel classification.  In _LastCFP: ICML Workshop on Learning from Multi-label data_ , 2010\.
  * Eyben et al. (2010) Eyben, F., Wöllmer, M., and Schuller, B.  Opensmile: the munich versatile and fast open-source audio feature extractor.  In _Proceedings of the 18th ACM international conference on Multimedia_ , pp. 1459–1462, 2010.
  * Eyben et al. (2015) Eyben, F., Scherer, K. R., Schuller, B. W., Sundberg, J., André, E., Busso, C., Devillers, L. Y., Epps, J., Laukka, P., Narayanan, S. S., et al.  The geneva minimalistic acoustic parameter set (gemaps) for voice research and affective computing.  _IEEE transactions on affective computing_ , 7(2):190–202, 2015.
  * Gales et al. (2014) Gales, M. J., Knill, K. M., Ragni, A., and Rath, S. P.  Speech recognition and keyword spotting for low-resource languages: Babel project research at cued.  In _Fourth International workshop on spoken language technologies for under-resourced languages (SLTU-2014)_ , pp. 16–23. International Speech Communication Association (ISCA), 2014.
  * Godfrey & Holliman (1997) Godfrey, J. J. and Holliman, E.  Switchboard-1 release 2.  _Linguistic Data Consortium, Philadelphia_ , 926:927, 1997\.
  * Hall et al. (2009) Hall, J. A., Andrzejewski, S. A., and Yopchick, J. E.  Psychosocial correlates of interpersonal sensitivity: A meta-analysis.  _Journal of nonverbal behavior_ , 33(3):149–180, 2009.
  * Hsu et al. (2021a) Hsu, W.-N., Bolte, B., Tsai, Y.-H. H., Lakhotia, K., Salakhutdinov, R., and Mohamed, A.  Hubert: Self-supervised speech representation learning by masked prediction of hidden units.  _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , 29:3451–3460, 2021a.
  * Hsu et al. (2021b) Hsu, W.-N., Sriram, A., Baevski, A., Likhomanenko, T., Xu, Q., Pratap, V., Kahn, J., Lee, A., Collobert, R., Synnaeve, G., et al.  Robust wav2vec 2.0: Analyzing domain shift in self-supervised pre-training.  _arXiv preprint arXiv:2104.01027_ , 2021b.
  * Ioffe & Szegedy (2015) Ioffe, S. and Szegedy, C.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In _International conference on machine learning_ , pp. 448–456. PMLR, 2015.
  * Kahn et al. (2020) Kahn, J., Rivière, M., Zheng, W., Kharitonov, E., Xu, Q., Mazaré, P. E., Karadayi, J., Liptchinsky, V., Collobert, R., Fuegen, C., Likhomanenko, T., Synnaeve, G., Joulin, A., Mohamed, A., and Dupoux, E.  Libri-light: A benchmark for asr with limited or no supervision.  In _ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pp. 7669–7673, 2020.  https://github.com/facebookresearch/libri-light.
  * Kharitonov et al. (2020) Kharitonov, E., Rivière, M., Synnaeve, G., Wolf, L., Mazaré, P.-E., Douze, M., and Dupoux, E.  Data augmenting contrastive learning of speech representations in the time domain.  _arXiv preprint arXiv:2007.00991_ , 2020.
  * Kingma & Ba (2015) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  In _Proc. ICLR_ , San Diego, USA, May 2015.
  * Lima et al. (2013) Lima, C. F., Castro, S. L., and Scott, S. K.  When voices get emotional: A corpus of nonverbal vocalizations for research on emotion processing.  _Behavior research methods_ , 45(4):1234–1245, 2013.
  * Maas et al. (2013) Maas, A. L., Hannun, A. Y., and Ng, A. Y.  Rectifier nonlinearities improve neural network acoustic models.  In _ICML Workshop on Deep Learning for Audio, Speech and Language Processing_ , 2013.
  * Nam et al. (2017) Nam, J., Loza Mencía, E., Kim, H. J., and Fürnkranz, J.  Maximizing subset accuracy with recurrent neural networks in multi-label classification.  _Advances in neural information processing systems_ , 30, 2017.
  * Panayotov et al. (2015) Panayotov, V., Chen, G., Povey, D., and Khudanpur, S.  Librispeech: an asr corpus based on public domain audio books.  In _2015 IEEE international conference on acoustics, speech and signal processing (ICASSP)_ , pp. 5206–5210. IEEE, 2015.
  * Pepino et al. (2021) Pepino, L., Riera, P., and Ferrer, L.  Emotion recognition from speech using wav2vec 2.0 embeddings.  _arXiv preprint arXiv:2104.03502_ , 2021.
  * Pratap et al. (2020) Pratap, V., Xu, Q., Sriram, A., Synnaeve, G., and Collobert, R.  Mls: A large-scale multilingual dataset for speech research.  _arXiv preprint arXiv:2012.03411_ , 2020.
  * Read et al. (2011) Read, J., Pfahringer, B., Holmes, G., and Frank, E.  Classifier chains for multi-label classification.  _Machine learning_ , 85(3):333–359, 2011.
  * Read et al. (2021) Read, J., Pfahringer, B., Holmes, G., and Frank, E.  Classifier chains: a review and perspectives.  _Journal of Artificial Intelligence Research_ , 70:683–718, 2021.
  * Saeki et al. (2022) Saeki, T., Xin, D., Nakata, W., Koriyama, T., Takamichi, S., and Saruwatari, H.  Utmos: Utokyo-sarulab system for voicemos challenge 2022.  _arXiv preprint arXiv:2204.02152_ , 2022.
  * Sauter et al. (2010) Sauter, D. A., Eisner, F., Ekman, P., and Scott, S. K.  Cross-cultural recognition of basic emotions through nonverbal emotional vocalizations.  _Proceedings of the National Academy of Sciences_ , 107(6):2408–2412, 2010.
  * Scherer (1994) Scherer, K. R.  Affect bursts.  _Emotions: Essays on emotion theory_ , 161:196, 1994.
  * Scherer & Scherer (2011) Scherer, K. R. and Scherer, U.  Assessing the ability to recognize facial and vocal expressions of emotion: Construction and validation of the emotion recognition index.  _Journal of Nonverbal Behavior_ , 35(4):305–326, 2011.
  * Schmitt & Schuller (2017) Schmitt, M. and Schuller, B.  Openxbow: introducing the passau open-source crossmodal bag-of-words toolkit.  _The Journal of Machine Learning Research_ , 18(1):3370–3374, 2017.
  * Schneider et al. (2019) Schneider, S., Baevski, A., Collobert, R., and Auli, M.  wav2vec: Unsupervised pre-training for speech recognition.  _arXiv preprint arXiv:1904.05862_ , 2019.
  * Schuller et al. (2016) Schuller, B., Steidl, S., Batliner, A., Hirschberg, J., Burgoon, J. K., Baird, A., Elkins, A., Zhang, Y., Coutinho, E., Evanini, K., et al.  The interspeech 2016 computational paralinguistics challenge: Deception, sincerity & native language.  In _17TH Annual Conference of the International Speech Communication Association (Interspeech 2016), Vols 1-5_ , pp. 2001–2005, 2016\.
  * Trouvain & Truong (2012) Trouvain, J. and Truong, K. P.  Comparing non-verbal vocalisations in conversational speech corpora.  In _Proceedings of the LREC Workshop on Corpora for Research on Emotion Sentiment and Social Signals_ , pp. 36–39. Citeseer, 2012.
  * Wagner et al. (2022) Wagner, J., Triantafyllopoulos, A., Wierstorf, H., Schmitt, M., Eyben, F., and Schuller, B. W.  Dawn of the transformer era in speech emotion recognition: closing the valence gap.  _arXiv preprint arXiv:2203.07378_ , 2022.
  * Wang et al. (2021) Wang, Y., Boumadane, A., and Heba, A.  A fine-tuned wav2vec 2.0/hubert benchmark for speech emotion recognition, speaker verification and spoken language understanding.  _arXiv preprint arXiv:2111.02735_ , 2021.


# Quasi-Random Frequency Sampling for Optical Turbulence Simulations

Amokrane Berdja1,2,* Massinissa Hadjara3,4  Marcel Carbillet5  Rafael L.
Bernardi6  and Romain G. Petrov5  1Centro de Astronomía (CITEVA), Universidad
de Antofagasta, Avenida Angamos 601, Antofagasta 1270300, Chile
2DICTUC S.A., Av. Vicuña Mackenna 4860, 7820436 Santiago, Chile
3Nanjing Institute of Astronomical Optics & Technology, National Astronomical
Observatories, Chinese Academy of Sciences, Nanjing 210042, China
4Chinese High Angular Resolution Southern Astronomical Laboratory
(CHARSAL)/Astro-Photonics Laboratory/Space and Planetary Exploration
Laboratory (SPEL), Dept. of Electrical Engineering, FCFM, Universidad de
Chile, Av. Tupper 2007, Santiago, Chile
5Laboratoire Lagrange (Université Côte d’Azur, Observatoire de la Côte d’Azur,
CNRS), Bât. Fizeau, Parc Valrose, 06100 Nice, France
6Department of Electrical Engineering, Pontificia Universidad Católica de
Chile, Vicuña Mackenna 4860, 7820436 Santiago, Chile
<EMAIL_ADDRESS>
††journal: opticajournal††articletype: Research Article

Optical turbulence modelling and simulation are crucial for developing
astronomical ground-based instruments, laser communication, laser metrology,
or any application where light propagates through a turbulent medium. In the
context of spectrum-based optical turbulence Monte-Carlo simulations, we
present an alternative approach to the methods based on the Fast Fourier
Transform (FFT) using a quasi-random frequency sampling heuristic. This
approach provides complete control over the spectral information expressed in
the simulated measurable, without the drawbacks encountered with FFT-based
methods such as high-frequency aliasing, low-frequency under-sampling, and
static sampling statistics. The methods’s heuristics, implementation, and an
application example from the study of differential piston fluctuations are
discussed.

##  1 Introduction

Optical turbulence numerical simulation has become an essential tool in
Astronomy, laser communication, and metrology. It is used to model how ground-
based instruments respond to optical turbulence in astronomy. For example, it
is central to the modelling of Adaptive Optics (AO) performance and
responsiveness [1, 2, 3]. It can also be used to test the effectiveness of any
instrument that estimates and monitors optical turbulence. Imaging
applications, such as in AO, rely on simulating phase fluctuations, including
scintillation, at a telescope’s pupil plane after propagation through the
atmosphere. Several simulation methods that have been proposed, including the
use of Zernike bases [4, 5, 6], but the most commonly used method is the
spectrum-based Monte-Carlo simulation. It relies on a model of the 2D spatial
power spectrum density of phase fluctuations (Kolmogorov, von Kármán,
Greenwood-Tarazano, "exponential", etc.), depending on how the effect of the
outer scale of optical turbulence is modelled [7]. The implementation
typically involves using the 2D Fast Fourier Transform (FFT) algorithm [8].

It has long been recognised that FFT-based simulations may suffer from certain
drawbacks in practice. This is mainly due to limited computational resources,
which can result in compromises in the sampling of spectra. The most well-
documented issue is the under-sampling of low frequencies in the case of phase
fluctuation simulations. Many approaches involving the addition of sub-
harmonics have been applied to address this issue. Another drawback is that
the spectral under-sampling is static. This means that the statistical
properties of the simulated quantity correspond to the averaged sampled power
spectrum density rather than the theoretical average power spectrum density
[9].

In addition to these drawbacks there is also an issue with many applications,
such as imaging through a circular aperture, where a significant portion of
the simulated phase fluctuations goes unused. This waste is even more
substantial in applications like stellar interferometry simulations when
phase-screens are used [10].

This paper proposes a solution for spectrum-based simulations of phase and
phase-related optical turbulence quantities. The solution relies on a quasi-
random sampling of the frequency space. This approach may solve many of the
inconveniences encountered with conventional FFT-based simulations when used
properly.

The following section introduces spectral methods and FFT-based simulations
(refer to Section 2). It also addresses issues that arise in practice when
using the FFT-based approach. We then discuss how our approach deals with
these issues (refer to Section 3). To demonstrate the statistical validity of
the method, we apply it to simulate the differential piston fluctuations
(refer to Section 4). This serves as a test for the statistical properties of
the simulated piston, as the structure function of these fluctuations can be
predicted analytically.

##  2 Spectrum-based and FFT-based simulations of optical turbulence

###  2.1 Spectrum-based simulation

The primary objective of of optical turbulence Monte-Carlo simulations is to
generate pseudo-random fields of quantities such as phase fluctuations to
replicate real measurements. The statistical properties of these fields are
determined by a given model.

The primary statistical properties of interest are second order statistics
such as the variance, the covariance, or the structure function. These
statistical properties are directly related to the power spectrum density of
the quantity being simulated as per the Wiener–Khinchin theorem. Therefore, a
fluctuating quantity $\chi$ with an average power spectrum density $W_{\chi}$,
which corresponds to a certain covariance/structure function, can be obtained
from a randomised spectrum whose squared modulus averages to $W_{\chi}$. The
following sections consider a 2D spatial field, which is commonly used in
optical turbulence simulations. However, the concepts discussed can be
generalised to temporal, 3D spatial, or spatial-temporal fields and power
spectrum densities. In this context, we are discussing the simulation of a
quantity $\chi(\vec{r})$, where $\vec{r}$ is a 2D spatial position vector, and
its corresponding theoretical spatial power spectrum density
$W_{\chi}(\vec{f})$, where $\vec{f}$ is a 2D spatial frequency vector. In this
case, each frequency $\vec{f}$ contributes a random harmonic to the field
$\chi(\vec{r})$.

Each random harmonic has a random amplitude and a random phase. As there are
no constraints on the phase of these random harmonics, it is appropriate to
assign them a uniform random distribution between $0$ and $2\pi$. In this
case, it can be stated that the phase as a function of spatial frequencies is
$2\pi R_{u}(\vec{f})$, where $R_{u}\sim\mathcal{U}(0,1)$ is a random (or
pseudo-random) number with a uniform distribution between $0$ and $1$.

Similarly, given that the only constraint on the random harmonic is to have a
finite variance, it is reasonable to postulate as per the maximum entropy
principle that the amplitude is normally distributed with a variance equal to
$W_{\chi}(\vec{f})$. Therefore, we can represent the amplitude of the random
harmonic as $R_{n}(\vec{f})\sqrt{W_{\chi}(\vec{f})}$, where
$R_{n}\sim\mathcal{N}(0,1)$ represents a random number with a normal
distribution and a variance of $1$. This assumption may need to be updated if
there are additional constrains on the spectrum from theory or observation.

With these considerations in mind, we can write a general encapsulation of
Monte Carlo spectral density-based simulation methods as:

$\chi(\vec{r})=\mathcal{F}^{-1}\left(R_{n}(\vec{f})\sqrt{W_{\chi}(\vec{f})}\exp{\left(2\pi
iR_{u}(\vec{f})\right)}\right).$ (1)

###  2.2 FFT-based implementation

Most implementations of the spectral-density approach rely on the numerical
efficiency of the FFT algorithm [8]. The 2D FFT implementations, especially
for phase fluctuations, are widely used in the simulation of isoplanatic and
anisoplanatic imaging applications [11, 12], as well as in AO in astronomy, to
name just a few examples. The use of 2D phase fluctuation matrices is very
convenient for introducing such effects as Fresnel diffraction [13] within the
Fourier optics framework [14] to simulate various light propagation through
turbulence phenomena such as scintillation [15].

Based on the above formulation, most 2D FFT implementations can be represented
as [16]:

$\chi(\vec{r})=\mathcal{FFT}^{-1}\left((R_{n}(\vec{f})+iR_{n}(\vec{f}))\sqrt{W_{\chi}(\vec{f})}\right).$
(2)

It is worth noting here that the two $R_{n}\sim\mathcal{N}(0,1)$ appearing in
Eq. 2 are independent.

This method can simulate not only phase fluctuations but also other optical
turbulence quantities such as angle of arrival fluctuations or the piston
fluctuations over a 2D spatial field, provided that the appropriate power
spectrum density is used. These quantities can often be expressed in the
Fourier space as a linear filtering of the phase fluctuations power spectrum
density. For instance, Fig. 1 displays the results of two 2D simulations: one
of piston fluctuations over a 2D field (right panel) and another of phase
fluctuations (left panel). Both simulations use the same atmospheric
turbulence parameters, and the impact of the telescope’s pupil filtering is
evident.

Figure 1: Phase (left) and piston (right) simulations on the same scale
obtained with the 2D FFT-based simulation technique.

###  2.3 Practical limitations

The FFT-based simulation can theoretically be applied to any case, as long as
the spatial frequency sampling is adjusted accordingly. However, due to
limited computational resources in terms of memory usage and execution time,
compromises are often made in frequency sampling. As a result, the spectral
information may not be optimally transferred to the simulated quantity.

The power spectrum densities of phase fluctuations are tipically a power-law
function in the form of $f^{-11/3}$. As a result, low frequencies are often
under-sampled, particularly at the origin $f=0$. This represents a frequency
domain of size $\delta\\!f\times\delta\\!f$ in the constant sampling required
for the FFT algorithm, with $\delta\\!f$ being the sampling pitch in the
frequency domain.

To address the issue of low-frequency problems [17], sub-harmonics are
incorporated to compensate for the information that is not captured by the 2D
FFT alone [18, 19, 20] and to break the FFT-induced circular symmetry. This
correction is adequate for some imaging applications, even though
incorporating sub-harmonics results in additional computational time penalties
[21]. Sub-harmonics can be easily implemented when dealing with phase
fluctuations spectra, but they become more cumbersome for other quantities
where sub-samplings do not necessarily occur at the origin of the frequency
space.

Overall, under-sampling, including sub-harmonics under-sampling, because of
the limited computational resources, leads to an approximation of the true
power spectrum density. This sampling, as it is static, introduces systematic
deviations from the desired theoretical ensemble statistics. As a consequence,
for a sequence of random phase (or other quantity) screens, the same set of
spatial frequencies are represented in each realisation, and the ensemble
statistical properties (covariance, structure function) correspond to those
sampled frequencies, not to the continuous spatial spectrum as defined by the
model [9].

One way to address the issue of static spatial frequencies is to randomise the
frequencies that contribute to the simulated quantity of each realisation.
This approach should result in ensemble statistics that correspond to a power
spectrum density covering the entire spectral continuum, as calculated from a
long sequence of realisations.

Several spatial frequency randomisation schemes have been proposed. One
approach is to randomise the frequencies completely [22], while another is to
apply a random jitter to a regular grid [23]. It is also possible to randomly
shift a regular grid [21], which can still be used with the 2D FFT in this
case. Additionally, randomised spatial frequencies can be introduced based on
the symmetry and distribution of the power spectrum density [24, 25, 9].

It is important to note that there is an additional consideration to take into
account, in addition to the constraints that define the statistics of every
spectral harmonic (refer to Section 2.1). To adhere to the energy cascade
model of fully developed turbulence, each simulation realisation should
ideally incorporate all spatial frequencies.

Spectral continuity is implicitly assumed in FFT-based simulations due to the
regular sampling, but this assumption may not hold true for random sampling.
To demonstrate this, one could sample a single spatial frequency for each
realisation, resulting in a single 2D harmonic with a random amplitude and
phase. The given realisation does not represent an optical turbulence
measurement, despite the ensemble statistics converging to the theoretical
model through averaging across multiple realisations. Ensemble statistics
alone are necessary but not sufficient to ensure accurate simulations.
Spectral continuity is also required.

A random or semi-random sampling of spatial frequencies could address the
statistical bias caused by the static sampling in FFT-based simulations and
introduce spectral continuity. However, it may also introduce additional
effects on the simulated fields due to discrepancy (refer to Section 3).

##  3 Low discrepancy quasi-random frequency sampling

In the following, Eq. 1 is simplified to the following form:

$\chi(\vec{r})=\frac{2}{\sqrt{2}}\sum_{f_{x},f_{y}}R_{n}(\vec{f})\sqrt{\delta\\!f_{x}\delta\\!f_{y}W_{\chi}(\vec{f})}\cos{\left(2\pi\left(\vec{f}\cdot\vec{r}+R_{u}(\vec{f})\right)\right)},$
(3)

where $f_{x}$ and $f_{y}$ are the orthogonal components of the 2D spatial
frequency vector $\vec{f}$.

Note that the factor of $2$ over the square root of $2$ in this equation is
due to the assumption that a sampled frequency vector $\vec{f}$ does not have
a counterpart vector $-\vec{f}$ as in the FFT method for instance. If a
counterpart vector were present, this factor would be $1$.

This approach does not take advantage of the optimised FFT routines, but it
offers the flexibility to simulate at specific areas of interest $\vec{r}$.
For example, if a simple model of piston fluctuations in an interferometric
system needs to be simulated, the effects are simulated precisely where they
are observed in the 2D space field. Similarly, for imaging applications, phase
fluctuations can be simulated inside an aperture or a collection of apertures,
with no need to simulate in the 2D space where light is masked. Another
advantage is the flexibility in sampling the spatial frequency components
required to simulate the desired quantity. For example, if the power spectrum
density hypothetically has a hard low-pass cutoff frequency, there is no need
to use frequencies beyond this limit. It is better to focus computational
effort on sampling parts of the spectrum that provide information about the
statistical properties of the simulated quantity.

###  3.1 Randomly jittered grid frequency sampling

As mentioned in Section 2.3, when sampling the spatial frequencies randomly
the ensemble statistical moments will inevitably converge, through a long
sequence of simulated sequences, to the theoretical model. If the frequencies
are randomly sampled across a spectral domain, a single simulated field may
have an uneven concentration of frequencies. The assumption of a fully
developed turbulence, upon which the theoretical models of optical turbulence
are built, implies however a uniform representation of harmonics throughout
the whole of the spectral domain.

To correct this imbalance in the spectral distribution, one solution would be
to maintain a regular grid $\vec{f^{\prime}}(f^{\prime}_{x},f^{\prime}_{y})$
with grid cells of size $\delta\\!f^{\prime}_{x}\times\delta\\!f^{\prime}_{y}$
to which a random jitter is introduced:

$\begin{cases}f_{x}=f^{\prime}_{x}+R_{u}(f_{x})\delta\\!f^{\prime}_{x},\\\
f_{y}=f^{\prime}_{y}+R_{u}(f_{y})\delta\\!f^{\prime}_{y}.\end{cases}$ (4)

In other words, if we subdivide the spectral domain into a grid of cells of
size $\delta\\!f^{\prime}_{x}\times\delta\\!f^{\prime}_{y}$, a spatial
frequency $\vec{f}(f_{x},f_{y})$ is randomly selected from within each cell.
$R_{u}\sim\mathcal{U}(0,1)$ are uniformly distributed.

This approach partially solves the large scale imbalance in the frequency
representation, but as with taking random frequencies from the whole domain,
there are many frequencies that will "clash", in the sense that they will be
too close together. It is as if the same frequency is used multiple times. The
left panel of Fig. 2 shows an example of a jittered regular grid where a few
spatial frequencies (red dots) are so close together that they are almost
indistinguishable.

When this happens, there are random gaps in the frequency coverage, but this
is not a matter of concern because the amplitude of each frequency is normally
distributed – this is naturally expected to happen anyway. However, in the
construction of the spectrum-based simulation approach, it is expected that
the probability of each harmonic component to be bounded by the normal
distribution that is assigned to it. If in the simulation of an optical
turbulence quantity, some frequencies have an excessively large amplitude,
even if the sequence converges statistically over many simulated fields, it is
not obvious whether the individual fields have the structure that mimics the
structure that can be observed in nature. Unless this is experimentally
observed or theoretically jusified, we want to avoid this so-called "clashing"
of frequencies.

Figure 2: Randomly jittered (left panel) and quasi-random (right panel)
frequency sampling for one realisation. $f_{x}$ and $f_{y}$ are the axes of
the 2D spatial frequency. The green lines represent a regular grid like in the
FFT-based simulations.

To avoid this clustering (or clashing) of spatial frequencies, we propose
instead to use low-discrepancy sequences (quasi-random numbers) to sample the
spatial frequencies.

###  3.2 Quasi-random (low-discrepancy) frequency sampling

We propose the use of low-discrepancy sequences, also known as quasi-random
sequences (in contrast to pseudo-random numbers) to sample spatial frequencies
in optical turbulence simulations as expressed in Eq. 3. Quasi-random
sequences are fully deterministic sequences that have the property of
progressively covering, as uniformly as possible, the entire domain on which
they are defined. Pseudo-random sequences are designed to give the impression
that a generated number has no "memory" of the numbers preceding it in the
sequence. In quasi-random sequences, however, each number avoids the preceding
number as much as possible. The net result is that pseudo-random sampling,
including its jittered grid variant, produces clusters and gaps in the spatial
domain. Quasi-random sampling, while deterministic, is not repetitive, and
covers the spectral range evenly without clusters (frequency clashes) and
gaps. They have exactly the desired characteristics outlined above.

To our knowledge, these sequences have mainly been used to overcome the so-
called "curse of dimensionality" in high-dimensional integration [26]. We
argue that they are a very useful tool for the numerical simulation of optical
turbulence phenomena, overcoming the difficulties encountered in previous
approaches.

Despite the many quasi-random sequences available in the literature and in the
scientific libraries (Halton, Hammersley, Sobol), we choose to experiment with
the $R_{2}$-sequence presented by Roberts [27] for two main reasons:

First, it is the simplest and easiest to implement in code. Second, as an
irrational fraction method it offers an attractive intuitive rational.

Consider a given periodic function $f_{1}(t)$ with period $T$ and a similar
function $f_{2}(t)$ with a period $\alpha T$. If the functions coincide at
time $t_{0}$, the next time they coincide will be at time
$t_{0}+nT=t_{0}+m\alpha T$ where $n$ and $m$ are integers and $\alpha=n/m$. If
$\alpha$ is irrational then the functions will never coincide again after
$t_{0}$. When the outputs are reduced (in the modulo sense) to a bounded
domain of length $T$, the output of $f_{2}(t)$ will continuously fill the
interval without repetition as $t$ becomes larger. Of course, this condition
is not sufficient to ensure low discrepancy, which for any such a sequence
needs to be tested or demonstrated.

The $R_{2}$-sequence can be expressed, to be consistent with the previous
formulations, as follows:

$\begin{cases}f_{x}=\delta\\!f_{x}n_{x}(s_{x}-0.5),\\\
f_{y}=\delta\\!f_{y}n_{y}(s_{y}-0.5),\end{cases}$ (5)

where $n_{x}$ and $n_{y}$ are natural numbers such that
$\delta\\!f_{x}n_{x}=2max(f_{x})$ and $\delta\\!f_{y}n_{y}=2max(f_{y})$, thus
defining the limits of the frequency domain.

Here, $s_{x}$ and $s_{y}$ are the 2D quasi-random numbers, between $0$ and
$1$, which in the context of the $R_{2}$-sequence approach are simultaneously
updated $n_{x}\times n_{y}$ times with the following additive recurrence:

$\begin{cases}s_{x}=\left\\{s_{x}+\frac{1}{\psi}\right\\},\\\
s_{y}=\left\\{s_{y}+\frac{1}{\psi^{2}}\right\\},\end{cases}$ (6)

where $\\{\\}$ denotes taking the fractional part. $\psi=1.324717957244746...$
is the so-called the plastic number (the positive root of the $x^{3}=x+1$
equation), an irrational number only limited by the machine’s floating-point
precision.

The right panel in Fig. 2 shows a realisation of a quasi-random sampling. It
is worth noting that contrarily to the randomly jittered realisation in the
panel above, there are no visible frequency clusters and gaps.

The initial seed for the quasi-random sequence can be set to
$(s_{x},s_{y})=(0,0)$ or $(s_{x},s_{y})=(0.5,0.5)$, according to Roberts [27],
but it could also be set to any random value according to our numerical
experiments. When simulating different realisations, it is important that the
quasi-random numbers are not reinitialised. One can draw an analogy with the
pseudo-random numbers that can be seeded once for the whole sequence of
simulations. In fact, $s_{x}$ and $s_{y}$ should be continuously updated
throughout the whole process, even when simulating a long series of fields,
because on the long run, the quasi-random frequencies will cover the whole of
the spectral range even so finely. It is precisely this property that allows a
faster convergence to more accurate statistical quantities.

In case the power spectrum density needs to be expressed in polar coordinates,
the quasi-random frequencies $(f_{x},f_{y})$ can be transformed into polar
quasi-random frequencies through a "square to disc" mapping [28].

Figure 3: Randomly jittered (left panel) and quasi-random (right panel)
frequency sampling for several realisations. Each colour represents a
different realisation. $f_{x}$ and $f_{y}$ are the axes of the 2D spatial
frequency.

Fig. 3 shows a comparison between, over many realisations of the simulation
(represented by different colours), the jittered grid sampling (left panel)
and the quasi-random sampling (right panel). The jitterd grid still has
clusters and gaps in the 2D frequency coverage. The quasi-random sampling, on
the other hand, in addition to being evenly distributed for each realisation
(each colour separately) covers the frequency space even more finely due to
its non-repeating nature, while maintaining the even distribution.

The formulation of the 2D quasi-random frequency sampling, as shown in Eq. 5,
can be simplified in most applications by taking
$\delta\\!f_{y}=\delta\\!f_{x}=\delta\\!f$ and $n_{x}=n_{y}=N$. This leaves us
with two free "hyper-parameters" $\delta\\!f$ and $N$. These two parameters,
under this formulation, give the flexibility to choose the spectral content of
the simulation. Once $\delta\\!f$ has been chosen, a value can be set for $N$,
where $max(f_{x,y})=N\delta\\!f$ being the highest frequency in either axis,
according to the power spectrum density and the required numerical precision
of the simulation. The next section presents an example of a simulation
describing this process.

##  4 Application Example

###  4.1 Piston fluctuations model

To demonstrate the practicality of using quasi-random sequences for optical
turbulence simulations, we present an example of simulating monochromatic
piston fluctuations for a long baseline stellar interferometer.

The simulation is designed to simulate a fringe-tracker [29, 30], which
measures the differential pistons in a specific wavelength band. This allows
for near-instantaneous corrections to be made in another band where the
"science object" is observed. The discussion at hand does not cover the
simulation details of the fringe tracker, which include polychromatic
dispersion, beam-combination, angular dispersion, anisoplanatism, measurement
errors, and fringe wander reconstruction.

The benefit of this option is its simplicity. If we were simulating imaging
through a telescope or a Michelson interferometer, we would need to discuss
minor technical details such as aperture sampling and vectorisation. However,
since this is not the case, we can exclude these details from our discussion.

In this example, there is a single turbulent layer that moves at a constant
velocity over the telescopes. According to the Taylor frozen turbulence
hypothesis, the temporal statistical properties are equivalent to the spatial
properties along the wind orientation. In other words, a time interval $\tau$
is equivalent to a distance of $x=\tau v_{w}$ where $v_{w}$ is the wind
velocity. $x$ is measured along the wind’s direction as represented in Fig. 5.
Although turbulence cannot remain frozen during a long time, for practical
purposes, we are only interested in the instrument’s response to differential
piston fluctuations (or equivalently, Optical Length Difference variations)
over times in the range of a few tens of milliseconds. Therefore, the Taylor
hypothesis is appropriate for this simulation.

A simplified model of the spatial power spectrum density of piston
fluctuations based on the von Kármán phase fluctuations model [7] is given by:

$W_{\overline{\varphi}}(\vec{f})\approx 0.02289\
r_{0}^{-\frac{5}{3}}\left(f^{2}+L_{0}^{-2}\right)^{-\frac{11}{6}}4\left(\frac{DJ_{1}(\pi
Df)-dJ_{1}(\pi df))}{\pi(D^{2}-d^{2})f}\right)^{2},$ (7)

where $r_{0}=r_{0}(\lambda)$ is the Fried parameter at wavelength $\lambda$,
$L_{0}$ is the outer scale of optical turbulence, $D$ is the telescope’s
diameter, and $d$ is the telescope’s central obstruction diameter.
$J_{1}(...)$ is the first kind Bessel function of order 1.

The von Kármán power spectrum density of phase fluctuations is represented by
the first term in Eq. 7 with the exponent $-11/6$. The equation’s second term,
which has an exponent of $2$, represents the aperture filtering. In this case,
the aperture is modelled as a circular aperture with a diameter of $D$ and a
circular central obstruction with a diameter of $d$.

Applying the procedure outlined in Eq. 3 and Eq. 5 is straightforward with
this model. However, the hyper-parameters $N=n_{x}=n_{y}$ and
$\delta\\!f=\delta\\!f_{x,y}$ require a choice to be made.

###  4.2 Hyper-parameters

The hyper-parameters $\delta\\!f$ and $N$ are selected based on the properties
of the power spectrum density $W_{\overline{\varphi}}(\vec{f})$ that need to
be expressed in the simulation.

Long-baseline interferometry is sensitive to piston fluctuations over
distances that are of the order of the outer scale $L_{0}$. To ensure
consistent sampling of low frequencies in the von Kármán’s model during
simulation, a reasonable value for $\delta\\!f$ is $\delta\\!f=0.25/L_{0}$.

To determine $N$, such that $2max(f_{x,y})=N\delta\\!f$, it is possible to
follow the following steps: Since in this case the power spectrum density is
isotropic ($W_{\overline{\varphi}}(\vec{f})=W_{\overline{\varphi}}(f)$), it is
possible to calculate (from the model) the quantity
$g(x)=\int_{0}^{x}df^{\prime}f^{\prime}W_{\overline{\varphi}}(f^{\prime})$.
Due to the finite variance of the von Kármán model and to the aperture
filtering of the telescope, the function $g(x)$ saturates and converges to a
finite variance $\sigma_{\overline{\varphi}}^{2}$ at
$f^{\prime}\rightarrow\infty$. A stopping criterion can be defined such that
$g(x)-g(x-\delta\\!f)<\epsilon$, where
$\epsilon<<\sigma_{\overline{\varphi}}^{2}$.

There are other possible criteria that depend on the required precision of the
simulated fluctuations. This provides additional flexibility to the simulation
approach, allowing for a choice between precision and speed.

###  4.3 Numerical Implementation

The implementation of the simulation algorithm is straightforward given a
numerical function `sqrt_w` that calculates the absolute square root of the
power spectrum density $W_{\overline{\varphi}}(\vec{f})$.

Fig. 4 is a pseudo-code of the simulation algorithm. It can easily be
implemented in any programming language.

Figure 4: A pseudo-code of the simulation procedure using quasi-random
numbers.

The power spectrum density function receives optical turbulence parameters
such as $r_{0}(\lambda)$ and $L_{0}$ as inputs. The vectors `xi` and `yi`
($i=1,2,...$) represent the 2D position components where pistons are
simulated.

The positions sampled depend on the simulation’s needs. In this case, they are
selected to achieve millisecond time sampling. Unlike in the FFT-based method,
this sampling is not bound to the sampling in the Fourier space, providing
greater modelling flexibility. It is evident from this pseudo-code that the
seeds `sx` and `sy` are global variables that are updated by the simulation
process. This guarantees that the entire spectral range, delimited here by
`df` and `N` is evenly sampled if the simulated pistons are to be used for
statistical characterisation.

Figure 5: A representation of the overall geometry for the simulation example
(left panel). The four blue circles are the four UT telescopes of the VLTI as
seen from above. The green lines represent the location of optical turbulence
as "seen" by the telescope through time. Their direction is the wind
direction. The light lines between the telescopes are the six interferometric
baselines. On the right panel a single realisation of the simulation generates
four simultaneous time series of piston fluctuations for each telescope.

The authors have implemented the code for these results in Python, and
accelerated it with the numba JIT compiler. The codes are available upon
request.

For the statistical tests presented in Section 4.4, the parameters used for
the simulation correspond to the VLTI telescopes as shown in Fig. 5 in the
left panel. VLTI is the interferometric mode of the European Southern
Observatory’s (ESO) Very Large Telescope (VLT) located at the Paranal
observatory in Chile. The VLTI combines, in the infrared wavelength ranges, 4
large VLT telescopes called UT (Unit Telescopes), which are fixed and have a
diameter of $D$ = 8m and $d$=1.1m central obstruction each, for a maximum
baseline of 130m. The VLTI also has an auxiliary array of 4 ATs (Auxiliary
Telescopes) with a diameter of $1.8m$ that can be moved to a large number of
stations on the Paranal platform with a maximum baseline of 200m. In the
simulation we are only concerned with piston fluctuations at the VLT
telescopes.

The right panel in Fig. 5 shows a single realisation of the simultaneous
simulated piston fluctuations at the 4 telescopes of the VLTI. The smoothness
of the curves is due to the aperture filtering of the telescopes and the wind
velocity as the turbulence layer moves above them. For the simulation in Fig.
5, as well as the one leading to the results in Fig. 6, wind speed is
$10.2m/s$ and the outer scale $L_{0}$=22m. The hyper-parameters $\delta
f=0.25/22m^{-1}$ and $f_{max}$ = 0.25$m^{-1}$, hence $N$=44. In other
conditions and with other optical turbulence quantities like phase
fluctuations, these parameters would be different.

###  4.4 Statistical test

The simulation’s geometry produces a single time series for every telescope,
enabling a straightforward statistical test. For each pair of telescopes $i$
and $j$, the structure function of the differential piston fluctuations
$D_{\Delta\overline{\varphi}}^{i,j}(\tau)$ is given by
$\left<\left|\Delta\overline{\varphi}_{i,j}(t)\Delta\overline{\varphi}_{i,j}(t+\tau)\right|^{2}\right>_{t}$.
$<...>_{t}$ represents an average over all the available values $t$ on the
same simulation realisation. $\Delta\overline{\varphi}_{i,j}(t)$ are the
differential piston fluctuations between telescopes $i$ and $j$, given by
$\overline{\varphi}_{i}(t)-\overline{\varphi}_{j}(t)$.

The length of the simulated piston fluctuations time series is finite.
Therefore, the statistical error, which depends on the number of points used
to for every separation $\tau$ and the coherence time, increases with $\tau$.
If there are $M$ points separated by $\delta\tau$,
$D_{\Delta\overline{\varphi}}^{i,j}(\delta\tau)$ is calculated using $M-1$
data pairs, resulting in lower statistical error, while and
$D_{\Delta\overline{\varphi}}^{i,j}((M-1)\delta\tau)$ is calculated only once,
resulting in a higher statistical error. If the one-realisation structure
function converges to the theoretical structure function for an infinite
length of the simulated time series, then the statistics can considered
"ergodic". This can be assured by using the spectral continuity constraint
with adequate sampling.

Structure functions from a large number of independent realisations are
averaged out. If the simulation statistics are good, the averaged structure
function will converge to the theoretical structure function. However, if the
simulation statistics are biased due to inadequate sampling, the averaged
structure function will converge to a function that differs from the
theoretical structure function.

In Fig. 5, the average structure function of the Optical Path Difference
(OPD), $D_{OPD}^{1,2}(\tau)$, defined as
$D_{OPD}^{i,j}(\tau)=\left(\frac{\lambda}{2\pi}\right)^{2}D_{\Delta\overline{\varphi}}^{i,j}(\tau)$,
is compared to the expected theoretical result. $\lambda$ in the wavelength of
the light. The average structure function has been obtained from 2000
independent realisations.

Figure 6: Average structure function of the Optical Path Difference (OPD) of
one of the six baselines (telescopes 1 and 2, separated by 56.5m) compared to
the corresponding theoretical structure function. The right panel is a detail
of the same result. With 2000 independent realisations, statistical errors are
visible at large intervals. With more realisations, the curves become
indistinguishable.

The structure function resulting from our simulation almost perfectly matches
the theoretical structure function. For large time delays $\tau$, there are
still some statistical errors (left panel of Fig. 6). These errors are
inevitable but they can be reduced with more realisations. At smaller $\tau$
values, the structure function is more precise (right panel in Fig. 6).

The primary objective of the simulation is to generate data that closely
resemble the expected measurements. The structure function confirms that the
information from the power spectrum density is accurately translated into the
simulated quantities. As discussed in Section 2.1, the structure function,
especially when calculated over a large number of realisation, is necessary to
assess the validity of the the simulation. However, the simulation’s validity
is not solely determined by the structure function. It is completed by the
physics-based constraint of spectral continuity. In this case, the model’s
properties are accurately captured.

##  5 Conclusion

We have presented an approach to simulate the observable fluctuations of
optical turbulence using a quasi-random sampling of spatial frequencies within
the framework of spectrum-based simulations.

We have also discussed how this approach allows us to overcome some of the
difficulties encountered with the more common FFT-based simulations. It could
be argued that it is unfortunate that we do not take advantage of the
optimisations offered by the FFT algorithm. However, depending on the
situation, this new approach not only does not suffer from the same
limitations as the FFT-based simulations do, but it also offers some
advantages.

Taking the example of imaging the output of a Michelson interferometer, the
FFT-based simulation method requires generating matrices that are large enough
to contain the baseline between the telescopes and sufficiently sample each
telescope. The phase fluctuations used for the image are only a small fraction
of the generated phase screen. The proposed quasi-random frequencies approach
generates phase fluctuations only where they are useful, at the telescopes’
apertures. Not using the FFT approach can be computationally advantageous in
this case.

We do not expect this approach to entirely existing methods. However, when
applicable, it offers clear advantages, provided that one knows how to adjust
the hyper-parameters for a given task. We have mentioned the possibility of
mapping the "square" quasi-random sampling to a polar quasi-random sampling
for special cases of expressing the power spectrum density. Systematising
hyper-parameters with a mapping that takes into account the properties of the
power spectrum density may be possible, but it is not necessary at this point.
Fixing these hyper-parameters on a case-by-case basis allows for more
flexibility.

The simulation approach’s validity has been demonstrated through differential
piston fluctuations. This method can be applied to other optical turbulence
measurements, such as phase or angle-of-arrival fluctuations, provided that
the power spectrum density of the measurement is known, either theoretically
or empirically.

Funding Do later. Acknowledgments Do later. Disclosures The authors declare no
conflicts of interest. Data availability Data underlying the results presented
in this paper are not publicly available at this time but may be obtained from
the authors upon reasonable request.

## References

  * [1] M. Carbillet, C. Vérinaud, B. Femenía, _et al._ , “Modelling astronomical adaptive optics - I. The software package caos,” Monthly Notices of the Royal Astronomical Society 356, 1263–1275 (2005).
  * [2] E. H. Por, S. Y. Haffert, V. M. Radhakrishnan, _et al._ , “High Contrast Imaging for Python (HCIPy): an open-source adaptive optics and coronagraph simulator,” in _Adaptive Optics Systems VI,_ vol. 10703 L. M. Close, L. Schreiber, and D. Schmidt, eds., International Society for Optics and Photonics (SPIE, 2018), p. 1070342.
  * [3] M. J. Townson, O. J. D. Farley, G. O. de Xivry, _et al._ , “Aotools: a python package for adaptive optics modelling and analysis,” Opt. Express 27, 31316–31329 (2019).
  * [4] N. A. Roddier, “Atmospheric wavefront simulation using Zernike polynomials,” Optical Engineering 29, 1174 – 1180 (1990).
  * [5] M. R. Whiteley, M. C. Roggemann, and B. M. Welsh, “Temporal properties of the zernike expansion coefficients of turbulence-induced phase aberrations for aperture and source motion,” J. Opt. Soc. Am. A 15, 993–1005 (1998).
  * [6] N. Chimitt and S. Chan, “Anisoplanatic optical turbulence simulation for near-continuous Cn2 profiles without wave propagation,” Optical Engineering 62, 078103 (2023).
  * [7] V. P. Lukin, “Outer scale of turbulence and its influence on fluctuations of optical waves,” Physics-Uspekhi 64, 280 (2021).
  * [8] R. N. Bracewell, _The Fourier transform and its applications_ (McGraw Hill, 2000), chap. 11.
  * [9] M. Charnotskii, “Comparison of four techniques for turbulent phase screens simulation,” J. Opt. Soc. Am. A 37, 738–747 (2020).
  * [10] I. L. Porro, T. Berkefeld, and C. Leinert, “Simulation of the effects of atmospheric turbulence on mid-infrared visibility measurements with the mid-infrared interferometric instrument for the very large telescope interferometer,” Appl. Opt. 39, 1643–1651 (2000).
  * [11] J. P. Bos and M. C. Roggemann, “Technique for simulating anisoplanatic image formation over long horizontal paths,” Optical Engineering 51, 101704 (2012).
  * [12] R. C. Hardie, M. A. Rucci, S. R. Bose-Pillai, _et al._ , “Modeling and simulation of multispectral imaging through anisoplanatic atmospheric optical turbulence,” Optical Engineering 61, 093102 (2022).
  * [13] J. D. Schmidt, J. A. Tellez, and G. J. Gbur, “Semi-analytic simulation of optical wave propagation through turbulence,” Appl. Opt. 61, 9439–9448 (2022).
  * [14] J. W. Goodman, _Introduction to Fourier optics_ (W. H. Freeman, 2005), chap. 4, McGraw-Hill physical and quantum electronics series.
  * [15] R. A. Johnston and R. G. Lane, “Modeling scintillation from an aperiodic kolmogorov phase screen,” Appl. Opt. 39, 4761–4769 (2000).
  * [16] J. D. Schmidt, _Numerical simulation of optical wave propagation with examples in MATLAB_ (SPIE, 2010), chap. 9.
  * [17] R. Frehlich, “Simulation of laser propagation in a turbulent atmosphere,” Appl. Opt. 39, 393–397 (2000).
  * [18] A. G. R G Lane and J. C. Dainty, “Simulation of a kolmogorov phase screen,” Waves in Random Media 2, 209–224 (1992).
  * [19] G. Sedmak, “Implementation of fast-fourier-transform-based simulations of extra-large atmospheric phase and scintillation screens,” Appl. Opt. 43, 4527–4538 (2004).
  * [20] M. Carbillet and A. Riccardi, “Numerical modeling of atmospherically perturbed phase screens: new solutions for classical fast fourier transform and zernike methods,” Appl. Opt. 49, G47–G52 (2010).
  * [21] D. A. Paulson, C. Wu, and C. C. Davis, “Randomized spectral sampling for efficient simulation of laser propagation through optical turbulence,” J. Opt. Soc. Am. B 36, 3249–3262 (2019).
  * [22] V. V. Voitsekhovich, D. Kouznetsov, V. G. Orlov, and S. Cuevas, “Method of random wave vectors in simulation of anisoplanatic effects,” Appl. Opt. 38, 3985–3992 (1999).
  * [23] A. Berdja, J. Borgnino, and A. Irbah, “Fresnel diffraction and polychromatic effects on angle-of-arrival fluctuations,” Journal of Optics A: Pure and Applied Optics 8, 244–251 (2006).
  * [24] M. Charnotskii, “Sparse spectrum model for a turbulent phase,” J. Opt. Soc. Am. A 30, 479–488 (2013).
  * [25] M. Charnotskii, “Statistics of the sparse spectrum turbulent phase,” J. Opt. Soc. Am. A 30, 2455–2465 (2013).
  * [26] H. Niederreiter, _Random Number Generation and Quasi-Monte Carlo Methods_ (Society for Industrial and Applied Mathematics, 1992), chap. 1, CBMS-NSF Regional Conference Series in Applied Mathematics.
  * [27] M. Roberts, “The unreasonable effectiveness of quasirandom sequences,” (2018). Last accessed 15 September 2023.
  * [28] P. Shirley and K. Chiu, “A low distortion map between disk and square,” Journal of Graphics Tools 2, 45–52 (1997).
  * [29] R. G. Petrov, F. Allouche, A. Boskri, _et al._ , “Hierarchical Fringe Tracking, sky coverage and AGNs at the VLTI,” in _Optical and Infrared Interferometry and Imaging VIII,_ vol. 12183 A. Mérand, S. Sallum, and J. Sanchez-Bermudez, eds., International Society for Optics and Photonics (SPIE, 2022), p. 121830I.
  * [30] M. Hadjara, R. G. Petrov, S. Lagarde, _et al._ , “VERMILION: Visitor Extension spectRal sub-Mid-Infrared Light InterferOmeter iNstrument,” in _Optical and Infrared Interferometry and Imaging VIII,_ vol. 12183 A. Mérand, S. Sallum, and J. Sanchez-Bermudez, eds., International Society for Optics and Photonics (SPIE, 2022), p. 121831V.