# A Software-Defined Channel Sounder for Industrial Environments with Fast
Time Variance

Niels Hendrik Fliedner, Dimitri Block, Uwe Meier  inIT - Institute of
Industrial Information Technologies
OWL University of Applied Sciences
Lemgo, Germany
Email: {niels.fliedner, dimitri.block<EMAIL_ADDRESS>
###### Abstract

Novel industrial wireless applications require wide-band, real-time channel
characterization due to complex multi-path propagation. Rapid machine motion
leads to fast time variance of the channel’s reflective behavior, which must
be captured for radio channel characterization. Additionally, inhomogeneous
radio channels demand highly flexible measurements. Existing approaches for
radio channel measurements either lack flexibility or wide-band, real-time
performance with fast time variance.

In this paper, we propose a correlative channel sounding approach utilizing a
software-defined architecture. The approach enables real-time, wide-band
measurements with fast time variance immune to active interference. The
desired performance is validated with a demanding industrial application
example.

ACI

    adjacent channel interference
CTI

    cross-technology interference
PDU

    protocol data unit
WCS

    wireless communication system
SWT

    standardized wireless technology
QoS

    quality of service
UWB

    ultra wide-band
VSG

    vector signal generator
RTSA

    real-time spectrum analyzer
ISI

    inter-symbol interference
PHY

    physical layer
MAC

    medium access control
DSSS

    direct-sequence spread spectrum
CSMA/CA

    carrier-sense multiple access with collision avoidance
OFDM

    orthogonal frequency-division multiplex
PSK

    phase-shift keying
WLAN

    wireless local area network
WSAN-FA

    wireless sensor and actuator network for factory automation
HART

    highway addressable remote transducer protocol
LR-WPAN

    low-rate wireless personal area network
BLE

    Bluetooth Low Energy
PNO

    Profibus Nutzerorganisation
IWLAN

    industrial WLAN
LTE

    3GPP long-term evolution
BA

    building automation
FA

    factory automation
PA

    process automation
IoT

    internet of things
ISM

    industrial, scientific and medical
ESD

    energy spectral density
PSD

    power spectral density
ESD

    energy spectral density
RF

    radio frequency
RSA

    real time spectrum analyzer
SNR

    signal-to-noise ratio
VSG

    vector signal generator
MR

    medium request
SIR

    signal-to-interference ratio
SINR

    signal-to-interference-plus-noise ratio
PAPR

    peak-to-average-power ratio
MIMO

    multiple-input and multiple-output
DC

    direct current
LO

    local oscillator
CW

    continuous-wave
FMCW

    frequency-modulated CW
PDP

    power delay profile
PLE

    path-loss exponent
AWGN

    additive white Gaussian noise
USRP

    universal software radio peripheral
SDR

    software-defined radio
OOT

    out-of-tree
DSP

    digital signal processing
HDR

    hardware-defined radio
PN

    pseudo-noise
MLS

    maximum length sequence
FZC

    Frank Zadoff Chu
PACF

    periodic auto-correlation function
ACF

    auto-correlation function
CCF

    cross-correlation function
PCCF

    periodic cross-correlation function
TOSM

    through-open-short-match
LOS

    line-of-sight
NLOS

    non-line-of-sight
OLOS

    obstructed line-of-sight
FFT

    fast Fourier transform
IQ

    in-phase and quadrature
STFT

    short-term Fourier transform
FPGA

    field-programmable gate array
LFSR

    linear-feedback shift register
AGC

    automatic gain control
GNSS

    global navigation satellite system
PMT

    polymorphic type
OS

    operating system
SBC

    single-board computer

##  I Introduction

The recent increase in performance requirements for devices in Industry 4.0
wireless applications demands the development of new wireless communication
systems [1]. Such developments require channel characterization in areas
previously not fully investigated [2].

Compared to office and outdoor environments, industrial environments behave
drastically worse. They are ranging from highly absorbent to highly reflective
environments [3]. Reflections cause multi-path propagation, which increase
frequency variance in factory automation (FA) environments [4]. Multi-path
propagation leads to radio channel time dispersion causing inter-symbol
interference (ISI) in WCS [2]. Characterization of channels with multi-path
propagation requires wide-band radio channel measurements.

Further, non-stationary multi-path channel conditions also impact FA
environments WCS [4]. Especially rapid machine motion causes fast time
variance of the channel’s reflective behavior. This challenges time-critical
closed-loop Industry 4.0 applications with latencies below 1 ms [1, 2]. In
order to address those fast changes, the radio channel measurement demands to
be fast.

Additionally, industrial environments are inhomogeneous, which ranges from
highly reflective, complex machine cells to uniform storage areas and
transportation paths. Such inhomogeneity leads to high spatial diversity
within FA environments [3]. Spatial diversity requires detailed spatial
channel characterization, and therefore dense spatial channel measurements.
Thus, for FA environments the measurement methodology has to be flexible,
scalable and ideally embedded in a low-cost hardware.

Few publications have shown implementation of real-time channel sounding with
correlative sequences. [5] use software-defined radios with National
Instruments’ LabVIEW to program the correlation function directly into the
field-programmable gate array (FPGA) for enhanced performance. Any changes to
the measurement setup have to be applied separate from the measurement via
FPGA programming, which limits the flexibility. The hardware-defined HIRATE
channel sounder introduced by [6] and utilized by [2] shows even better
performance, but conversely is more limited to changes during measurement than
[5]. To the best of our knowledge, so far no publication implements a
software-defined real-time architecture for industrial radio channel
measurements.

Figure 1: Flexible and scalable architecture of the software-defined wireless
channel sounding approach corrsounder for industrial real-time channel
characterization

In this paper we propose a novel software-defined wireless channel sounding
approach named corrsounder for industrial real-time channel characterization.
In general, channel sounding techniques enable wide-band channel measurements,
and therefore enable multi-path propagation characterization. While channel
sounding suffers from limited accuracy, its periodic sequence-based
stimulation enable the characterization of fast varying environments.
Additionally, the software-defined realization enables flexible and scalable
implementation with off-the-shelf radio frequency (RF) front ends as shown in
Fig. 1. Furthermore, the proposed real-time capability allows for embedding
the measurement system into live processing on demand.

This paper is structured as follows: Section II provides an overview of
current solutions for industrial wireless channel sounding approaches, Section
III describes the the mathematical basics of the utilized correlative channel
sounding, Section IV details the proposed software architecture of
corrsounder, Section V validates the desired performance with an application
example, and finally a conclusion is drawn in Section VI.

##  II Industrial Wireless Channel Sounding

Figure 2: Signal and control flow graph of the sequence-based channel sounding
approach

For channel measurements there exist three methods of channel stimulation:
continuous-wave (CW), impulse-based and sequence-based stimulation. But CW\-
and impulse-based channel characterizations depend on emission measurements.
Hence, they are sensitive to interference. Therefore, using these methods it
is not possible to measure only the frequency or impulse response of an active
industrial environment.

Sequence-based stimulations overcome the time limitation and take advantage of
the repetitive structure. Thereby, various sequences exist for stimulation.
For reception, sequence-based stimulation require correlation methods, which
directly return the scaled impulse response. These correlation methods
suppress the impact of any active interference dissimilar to the utilized
sequence. In contrast to CW\- and impulse-based stimulations, these methods
enable the segregated measurement of the frequency or impulse response of an
industrial environment with additional transmitters, which is called an active
environment.

Generally speaking, correlative sequences can be differentiated into binary
and non-binary ones. Typical binary sequences are pseudo-noise (PN) sequences.
In general, binary sequences suffer from an inconstant spectral envelope, and
therefore linear distortion. They are commonly realized in the form of maximum
length sequences, which offer the highest sequence length $N_{\text{seq}}$ at
a certain linear-feedback shift register (LFSR) length $l$, with
$N_{\text{seq}}=2^{l}-1$. The periodic auto-correlation function (PACF) of
MLSs have a non-zero out-of-phase value $v_{\text{oop}}=-1$ [9], and thus add
a systematic error to the correlation, which limits the dynamic range to
$D_{\text{r}}=N_{\text{seq}}-|v_{\text{oop}}|$ [9]. Hence, they are not
considered perfect sequences. MLSs provide a constant magnitude, and hence a
peak-to-average-power ratio (PAPR) of one. This lowers the requirements on the
automatic gain control (AGC) of RF frontends, and therefore SDRs.

For non-binary sequences, there are perfect ones which provide the desired
property of a zero out-of-phase PACF value and therefore leads to the full
dynamic range of $D_{\text{r}}=N_{\text{seq}}$. An optimal choice for a non-
binary perfect sequence is called Frank Zadoff Chu (FZC) sequence [10]. They
also provide a constant magnitude, and hence a PAPR of one [11]. Additionally,
the spectral envelope of FZC sequences is constant without linear distortion.
Furthermore, these sequences are invulnerable to carrier frequency offset
[11], which allows more flexible radio channel measurements in complex
scenarios, since no further frequency synchronizing is needed. Additionally,
multiple FZC same-length sequences can be generated, which enables multiple-
input and multiple-output (MIMO) channel measurement [12] and multi-user time
synchronization, e.g. in 3GPP long-term evolution (LTE) communication [13].

##  III System Description

The sequence-based sounding signal flow is illustrated in Fig. 2.

For the time-discrete channel stimulation $x[t]$, a signal source outputs
infinite repetitions of a chosen correlative sequence $x_{\text{seq}}[t]$:

$\displaystyle x[t]$
$\displaystyle=\sum\limits_{i=-\infty}^{\infty}x_{\text{seq}}[t-i\,T_{\text{seq}}]$
(1)

with the duration $T_{\text{seq}}=T_{\text{s}}\cdot N_{\text{seq}}$, with the
sampling time $T_{\text{s}}$ and the length $N_{\text{seq}}$ . Then, the
resulting repetitive signal is transmitted with a SDR.

The transmitting and receiving SDRs involve in-phase and quadrature (IQ)
modulation and demodulation, respectivly. Additionally, they perform frequency
translation. Hence, the SDRs handle the low-pass equivalent of a band-pass
signal. So, a band-pass radio channel is stimulated and therfore also
evaluated.

The radio channel influences the transmitted signal:

$\displaystyle y[t]$ $\displaystyle=x[t]~{}\textasteriskcentered~{}h[t]+n$ (2)

with the equivalent low-pass impulse response $h[t]$ and additive white
Gaussian noise (AWGN) $n$. It is important to note, that $h[t]$ contains the
impact of the radio channel $h_{\text{ch}}[t]$ and also the linear distortion
caused by utilized antenna equipment such as cables $h_{\text{cable}}[t]$:

$\displaystyle h[t]$
$\displaystyle=h_{\text{ch}}[t]~{}\textasteriskcentered~{}h_{\text{cable}}[t]$
(3)

The following signal processing block, called sequence gate, ensures the
repetitive structure of the received signal $y[t]$ for correlation later-on.
The repetitions may be disrupted by trigger events such as reception buffer
overflow or from other external sources e.g. global navigation satellite
system (GNSS) localization events. The sequence gate mutes disrupted sequence
durations and passes only full ones. Its output signal can be expressed as:

$y_{\text{g}}[t]=\sum\limits_{j=-\infty}^{\infty}y[t-k_{j}\,T_{\text{seq}}]\cdot
w_{\text{rect}}\left({(t-k_{j}\,T_{\text{seq}})}/{T_{\text{seq}}}\right)$ (4)

with the strictly monotonous growing index $k_{j}>k_{j-1}$ and the rectangular
windowing function $w_{\text{rect}}(t/T)$, which is one for the time interval
$\left[0,T\right)$ and otherwise zero. With no trigger events,
$y_{\text{g}}[t]$ is equal to $y[t]$.

Then, the gated signal $y_{\text{g}}[t]$ is correlated with the chosen
stimulation sequence $x_{\text{seq}}[t]$, which can be expressed as conjugate
complex time-reversed convolution:

$\displaystyle y_{\text{corr}}[t]$
$\displaystyle=y_{\text{g}}[t]~{}\textasteriskcentered~{}x_{\text{seq}}^{*}[-(t-T_{\text{seq}})]$
(5)

With no trigger events, $y_{\text{corr}}[t]$ can be separated into a noise-
impact term $y_{\text{corr,n}}[t]$ and a signal term $y_{\text{corr,s}}[t]$:

$\displaystyle y_{\text{corr}}[t]=\,$
$\displaystyle(x[t]~{}\textasteriskcentered~{}h[t])~{}\textasteriskcentered~{}x_{\text{seq}}^{*}[-(t-T_{\text{seq}})]$
$\displaystyle+n~{}\textasteriskcentered~{}x_{\text{seq}}^{*}[-(t-T_{\text{seq}})]$
(6) $\displaystyle=\,$ $\displaystyle
y_{\text{corr,s}}[t]+y_{\text{corr,n}}[t]$ (7)

With good AWGN cross-correlation properties of the chosen stimulation sequence
as discussed in Section II, the impact of $y_{\text{corr,n}}[t]$ can be
neglected. Then, the associative property of convolution is applied:

$\displaystyle y_{\text{corr}}[t]$
$\displaystyle=(\,\underbrace{x[t]~{}\textasteriskcentered~{}x_{\text{seq}}^{*}[-(t-T_{\text{seq}})]}_{\text{\acs{PACF}}}\,)~{}\textasteriskcentered~{}h[t]$
(8)

The PACF property of correlative sequences with constant out-of-phase PACF
value $v_{\text{oop}}$ results in:

$\displaystyle y_{\text{corr}}[t]$
$\displaystyle=\Big{(}v_{\text{oop}}+(N_{\text{seq}}-v_{\text{oop}})\sum\limits_{i=-\infty}^{\infty}\delta[t-t_{i}]\Big{)}\textasteriskcentered
h[t]$ (9)

with the measurement time substitution
$t_{i}=(i+1)T_{\text{seq}}-T_{\text{s}}$.

Next, the distributive property of convolution and the sifting property of the
Dirac function are applied. The result is a time-shifted and value-scaled
impulse response series.

$\displaystyle y_{\text{corr}}[t]=\,$
$\displaystyle(N_{\text{seq}}-v_{\text{oop}})\sum\limits_{i=-\infty}^{\infty}h[t-t_{i}]$
(10) $\displaystyle+v_{\text{oop}}\,\textasteriskcentered\,h[t]$ (11)

Hence, a non-zero $v_{\text{oop}}$ adds a bias error
$v_{\text{oop}}\,\textasteriskcentered\,h[t]$, which is the case for non-
perfect sequences.

In order to return the desired impulse response $h[t]$, the signal
$y_{\text{corr}}[t]$ needs to be scaled according to the number of samples per
sequence $N_{\text{seq}}$.

$\displaystyle y_{\text{norm}}[t]$
$\displaystyle=\frac{1}{N_{\text{seq}}}\cdot y_{\text{corr}}[t]$ (12)

This also shows, that the bias error
$v_{\text{oop}}\,\textasteriskcentered\,h[t]$ gets reduced with increased
sequence length $N_{\text{seq}}$.

The forward transmission tracking error terms are used to correct the linear
distortion caused by $h_{\text{cable}}[t]$, which is still contained within
$y_{\text{norm}}[t]$. This correction is executed by convolution with
$h_{\text{ftt}}[t]$.

$\displaystyle y_{\text{ch}}[t]$
$\displaystyle=y_{\text{norm}}[t]~{}\textasteriskcentered~{}h_{\text{ftt}}[t]$
(13)
$\displaystyle\approx\sum\limits_{i=-\infty}^{\infty}h_{\text{ch}}[t-t_{i}]=\sum\limits_{i=-\infty}^{\infty}h_{\text{ch}}[\tau,t_{i}]$
(14)

with the finely granulated impulse response time $\tau$.

Thus, the proposed measurement method returns the sought-for channel impulse
response series.

##  IV Software Architecture

The channel sounding approach corrsounder is realized based on a modular
software architecture. In the following sub-sections the software architecture
shown in Fig. 1 is discussed.

###  IV-A Distributed Processing

Channel sounding with sequence-based stimulations can be splitted into the
transmission part for stimulation and the reception part for correlation.
These are independent from each other as shown in Fig. 2, as long as an
identical sequence is provided. Hence, channel sounding stimulation and
correlation may be separated, and the associated SDRs can be connected to
multiple processing units such as embedded host computers. Each processing
unit has to be connected to at least one SDR, performing either the
stimulation or correlation processing. With multiple SDRs it may perform any
arbitrary combination of the stimulation or correlation processing. Therefore,
the software architecture approaches a distributed processing. It also enables
integration into existing utilized wireless infrastructure, and therefore
increases flexibility.

###  IV-B Operating System

The corrsounder expects a Linux-based operating systems for the processing
units. Hence, general-purpose Ubuntu OS for a high-power x86 hardware
architecture, but also low-power embedded ARM-based Linux derivates for
single-board computer (SBC) can be addressed. The OSs have to provide real-
time hardware bindings to the SDRs such as device drivers. In case of
additional trigger sources, the bindings such as hardware clocks and GNSS
receivers also have to be provided by the OS . Additionally, the OS has to
provide the Python runtime and the computation library numpy for signal
processing.

###  IV-C Kernel Module

The kernel module pycorrsounder contains basic routines for channel sounding
operation as well as for evaluation. It only requires the Python related
libraries. Thus, the dependencies of corrsounder for pre- and post-processing
are limited up to the kernel module, which increases the flexibility.
pycorrsounder contains the following routines:

#### Sequence Generation and Correlation

Generation routines are provided for binary MLS and non-binary FZC sequences.
Additional sequence generation routines can be added later on. The module also
contains routines for sequence correlation: auto-correlation function (ACF),
PACF, cross-correlation function (CCF), and periodic cross-correlation
function (PCCF). The ACF can be used for evaluation of a finite signal, while
PACF is the counterpart for assumed infinite repetition of the finite signal.
The CCF and PCCF can be used accordingly for evaluation of two non-equal
signals, such as a transmitted sequence and the corresponding received channel
response.

#### Error Correction and Filtering

For two-port systematic linear error correction with two SDR antenna ports, a
through adjustment routine is provided. It requires the error term for forward
transmission tracking. The error term can be acquired by directly connecting
the SDR antenna ports. Through adjustment is easy to perform, removes the main
frequency response error and can be utilized when medium accuracy is required.
Various other distortion effects can be suppressed by the provided discrete
filtering routines. A removal routine for peaks and a down-sampling routine
are provided. The removal routine can suppress spectral and temporal peaks
such as the direct current (DC) bias. The bias is caused by the local
oscillator (LO) RF leakage of a receiving SDR [14]. It results in a constant
value offset and therefore a spectral peak around 0 Hz. Hence, the removal
routine fades out such spectral frequency range and interpolates the values
afterwards. The down-sampling routine is utilized to suppress high-frequency
attenuation by extracting the low-frequency part. The routine determines a
cutoff frequency at which the power spectral density (PSD) falls below a
threshold relative to the maximum, e.g. -3 dB for the edge frequency.

#### Channel Characterization

Routines are provided for channel characterization. Thereby, impulse response
time $\tau$ and measurement time $t$ evaluation routines are provided for
single-sequence and inter-sequence characterization, respectively. The former
characterize the frequency and the impulse response at a certain measurement
time $t$. Thereby, the power delay profile (PDP) $E[|h(\tau)|^{2}]$, the mean
delay $\bar{\tau}$ and the delay spread $\tau_{\text{rms}}$ can be
characterized. For the frequency response, the routines characterize the PSD
$E[|H(f)|^{2}]$ and statistical metrics like percentile 10 % $H_{10}$, 90 %
$H_{90}$, 50 % $H_{\text{median}}$ of the magnitude $|H(f)|$ and the coherence
bandwidth $B_{C}$. Additionally, arbitrary distributions can be estimated like
the Rician distribution for line-of-sight (LOS) conditions. Furthermore, the
path-loss exponent (PLE) can be estimated based on the transmission and
reception antenna gain. Measurement time characterization routines include for
example the Doppler spread $h[\tau,f_{d}]$ and the coherence time $T_{C}$.

###  IV-D GNU Radio

In case the processing unit has to perform real-time processing, the well-
known SDR framework GNU Radio is required. GNU Radio is a library of digital
signal processing (DSP) blocks, SDR bindings and application routines with
several command-line and graphical user interfaces. It is utilized in several
simulation studies as well as for DSP application deployment. For the signal
and control flow Fig. 2 different GNU Radio approaches are applied. For the
signal flow, stream-based items of complex IQ samples have to be utilized. GNU
Radio provides in-tree DSP blocks for stream processing such as SDR source and
sink, convolutional and DFT-based filter blocks. In contrast, the control flow
can also be event-based for triggers. Therefore, GNU Radio provides in-tree
blocks for event-based processing with polymorphic type (PMT) tags and PMT
messages. PMT tags are synchronous events which annotate signal streams. PMT
messages are asynchronous events which are processed concurrently and
independent of the signal streams.

###  IV-E corrsounder DSP Blocks

The corrsounder DSP blocks are provided within the out-of-tree (OOT) GNU Radio
module gr-corrsounder. It requires a GNU Radio instance and the routines of
the kernel module pycorrsounder. gr-corrsounder includes the DSP blocks
sequence source, sequence gate, sequence correlation, error correction,
channel characterization, and the meta information. These blocks are already
discussed in Section III.

###  IV-F corrsounder DSP Applications

The module gr-corrsounder also contains the corrsounder DSP applications.
These are executable flow graphs of connected DSP blocks. Mainly, four
applications are provided. The first application is called Channel Sounding.
It represents the channel sounding signal and control flow shown in Fig. 2. It
combines two applications: Sequence Stimulation and Sequence Correlation.
Sequence Stimulation is an application for the stimulation part with the
transmission-related processing blocks. Sequence Stimulation is the
counterpart for reception-related part. The final application is called
Calibration. It is a derivate of the Channel Sounding application without the
error correction blocks. It has to be utilized for calibration to determine
the error terms.

###  IV-G Source Code Availability

The source code of the channel sounder corrsounder is online available via the
URL https://github.com/inIT-HF/gr-corrsounder. It contains the module gr-
corrsounder with the DSP blocks and applications as well as the kernel module
pycorrsounder.

##  V Application Example

An exemplary measurement demonstrates the performance of our proposed channel
sounder architecture. The radio channel consists of a robot arm performing a
pick-and-place routine inside a separated cell. Two antennas are placed in
opposite corners with a distance of 3.8 m and a height of 1.5 m. The robot arm
obstructs the LOS periodically during motion.

We utilize distributed processing, with a vector signal generator (VSG) for
stimulation and a host computer connected to an SDR for correlation. The
stimulation sequence is generated by the kernel module pycorrsounder. For
correlation, the corrsounder DSP application Sequence Correlation is used.

TABLE I: Measurement Setup and Parameters Unit | Value
---|---
VSG | Rohde& Schwarz SMBV100A
Antennas | WiMo 18720.3H 5 GHz
Antenna gain | 5 dBi
Robot arm | Reis RV16
SDR | Ettus Research USRP X300
Host computer OS | Ubuntu Linux 14.04
Center frequency | 5.80 GHz
Sampling rate | 100 MSps
Sequence type | FZC
Sequence length | 1024 samples
FZC parameter | 7
Measurement 3 dB bandwidth | 93.07 MHz
Calibration routine | Through
DC bias suppression bandwidth | 781 kHz
Measurement duration | 12 s

With the utilized hardware listed in Table I, the proposed channel sounder
architecture enables real-time measurement without buffer overflows with the
mentioned performance.

The motion of the robot arm causes a frequency response pattern with a
repetition time of approximately 7.5 s, as shown in Fig. 3. In contrast, the
measurement time resolution of 10.24 µs falls far below the time variance of
the fast motion. Therefore, our proposed approach surpasses the fast time
variance requirements for industrial applications with even higher demands.

Figure 3: Frequency response

The Doppler plot in Fig. 4 shows the impulse response and indicates the time
variance of the multi-path propagation. The plot resolves distinct multi-path
components in the finely granulated impulse response time resolution of 10.75
ns. Furthermore, the long-term measurement enables Doppler spread resolution
of 0.083 Hz for time-sensitive characterization of the multi-path propagation.
Additionally, the maximum resolvable Doppler spread of 48.8 kHz corresponds to
a speed of 2526 m/s.

Figure 4: Doppler plot

The application example measurement results in a dynamic range of 45.1 dB.
This dynamic range enables a measurement distance of 695.87 m under ideal
propagation conditions.

The application example validates, that the proposed channel sounder
architectures is suited for wide-band, real-time radio channel measurements
with fast time variance.

##  VI Conclusion

Novel industrial wireless applications require wide-band, real-time channel
characterization due to complex multi-path propagation. Fast time variance of
the channel’s reflective behavior must be captured for radio channel
characterization. Additionally, inhomogeneous radio channels demand highly
flexible measurement measurements. Existing approaches for radio channel
measurements either lack flexibility or wide-band, real-time performance with
fast time variance.

The proposed approach utilizes a software-defined architecture for flexibility
and scalability. Furthermore, the approach uses correlative channel sounding
for real-time, wide-band measurements immune to active interference.
Additionally, the approach enables distributed processing, and therefore it is
favorable for embedded implementation.

In an application example, the distributed processing property is validated
for characterization of a radio channel obstructed by robot arm motion. Hence,
the proposed correlative channel sounding approach is suited for industrial
wide-band, real-time channel characterization with fast time variance.

## Acknowledgment

Parts of this work are funded by the German BMBF (project HiFlecs, 16KIS0266).

## References

  * [1] A. Frotzscher, U. Wetzker, M. Bauer, M. Rentschler, M. Beyer, S. Elspass, and H. Klessig, “Requirements and current solutions of wireless communication in industrial automation,” in _2014 ICC - 2014 IEEE International Conference on Communication Workshop (ICC)_ , 2014, pp. 67–72.
  * [2] B. Holfeld, D. Wieruch, L. Raschkowski, and T. Wirth, “Radio Channel Characterization at 5.85 GHz for Wireless M2M Communication of Industrial Robots,” in _WCNC 2016 - IEEE Wireless Communications and Networking Conference_ , 2016.
  * [3] P. Stenumgaard, J. Chilo, J. Ferrer-Coll, and P. Angskog, “Challenges and conditions for wireless machine-to-machine communications in industrial environments,” _IEEE Communications Magazine_ , vol. 51, no. 6, pp. 187–192, 2013.
  * [4] L. Tang, K.-C. Wang, and Y. Huang, “Performance Evaluation and Reliable Implementation of Data Transmission for Wireless Sensors on Rotating Mechanical Structures,” _Structural Health Monitoring_ , vol. 8, no. 2, pp. 113–124, 2009.
  * [5] J. Elofsson and P. Seimar, “Software Defined Radio based MIMO channel sounding,” Master’s Thesis, Lund University, Lund, May 25, 2016. [Online]. Available: http://lup.lub.lu.se/luur/download?func=downloadFile&recordOId=8874710&fileOId=8874711
  * [6] W. Keusgen, A. Kortke, M. Peter, and R. Weiler, “A highly flexible digital radio testbed and 60 GHz application examples,” 2013.
  * [7] Y. Konishi, M. Kim, M. Ghoraishi, J. Takada, S. Suyama, and H. Suzuki, “Channel Sounding Technique using MIMO Software Radio Architecture,” in _Proceedings of the 5th European Conference on Antennas and Propagation (EUCAP)_. Piscataway, NJ: IEEE, 2011, pp. 2546–2550.
  * [8] S. Boyd, “Multitone signals with low crest factor,” _IEEE Transactions on Circuits and Systems_ , vol. 33, no. 10, pp. 1018–1022, 1986\.
  * [9] G. Martin, “Wideband channel sounding dynamic range using a sliding correlator,” in _VTC 2000 - Spring_. Piscataway, NJ: IEEE, 2000, pp. 2517–2521.
  * [10] D. Sarwate, “Bounds on crosscorrelation and autocorrelation of sequences (Corresp.),” _IEEE Transactions on Information Theory_ , vol. 25, no. 6, pp. 720–724, 1979.
  * [11] H. Cao, Z. Zhao, W. Ni, M. El-Hadidy, and T. Kaiser, “Measurement and ray-tracing of wideband indoor channel in UHF TV white space,” in _2011 CogART - Proceedings of the 4th International Conference on Cognitive Radio and Advanced Spectrum Management_. New York, NY: ACM, 2011, pp. 1–5.
  * [12] J.-K. Hwang, Y.-T. Tsai, and J.-D. Li, “Multiuser channel estimation for SC-FDMA system with CAZAC sounding reference sequence and Instruments-in-MATLAB verification,” in _ISCE 2012 - IEEE 16th International Symposium on Consumer Electronics_. [Piscataway, NJ]: [IEEE], 2012, pp. 1–5.
  * [13] ETSI, “LTE; Evolved Universal Terrestrial Radio Access (E-UTRA); Physical channels and modulation (3GPP TS 36.211 version 12.5.0 Release 12),” 2015.
  * [14] R. Svitek and S. Raman, “DC offsets in direct-conversion receivers: Characterization and implications,” _IEEE Microwave Magazine_ , vol. 6, no. 3, pp. 76–86, 2005.

  *[FA]: factory automation
  *[ISI]: inter-symbol interference
  *[WCS]: wireless communication system
  *[FPGA]: field-programmable gate array
  *[RF]: radio frequency
  *[CW]: continuous-wave
  *[PN]: pseudo-noise
  *[LFSR]: linear-feedback shift register
  *[PACF]: periodic auto-correlation function
  *[MLSs]: maximum length sequence
  *[PAPR]: peak-to-average-power ratio
  *[AGC]: automatic gain control
  *[SDRs]: software-defined radio
  *[FZC]: Frank Zadoff Chu
  *[MIMO]: multiple-input and multiple-output
  *[LTE]: 3GPP long-term evolution
  *[SDR]: software-defined radio
  *[IQ]: in-phase and quadrature
  *[AWGN]: additive white Gaussian noise
  *[GNSS]: global navigation satellite system
  *[OS]: operating system
  *[SBC]: single-board computer
  *[OSs]: operating system
  *[MLS]: maximum length sequence
  *[ACF]: auto-correlation function
  *[CCF]: cross-correlation function
  *[PCCF]: periodic cross-correlation function
  *[DC]: direct current
  *[LO]: local oscillator
  *[PSD]: power spectral density
  *[PDP]: power delay profile
  *[LOS]: line-of-sight
  *[PLE]: path-loss exponent
  *[DSP]: digital signal processing
  *[PMT]: polymorphic type
  *[OOT]: out-of-tree
  *[VSG]: vector signal generator


$\boldsymbol{u}=(u_{1},u_{2},u_{3})^{T}\in[0,1]^{3}$, consider the copulas

$C_{\theta}(\boldsymbol{u})=u_{1}u_{2}u_{3}\\{1+\theta(1-u_{1})(1-u_{2})\\};$
(I.1a)
$\tilde{C}_{\theta}(\boldsymbol{u})=u_{1}u_{2}u_{3}\\{1+\theta(1-u_{1})(1-u_{2})(1-u_{3})\\},$
(I.1b)

where $|\theta|\leq 1$ (see Johnson and Kotz, 1975; Georges et al., 2001,
Remark 4, p. 9). Obviously, both $C_{\theta=0}$ and $\tilde{C}_{\theta=0}$
correspond to the independent copula case $u_{1}u_{2}u_{3}$.

Suppose $\boldsymbol{U}=(U_{1},U_{2},U_{3})^{T}\sim C_{\theta}$, see (I.1a).
Obviously, $(U_{1},U_{2})^{T}$ and $U_{3}$ are independent; also, one can
easily verify that the cdf of $U_{1}|U_{2}=u_{2}$ is
$F_{U_{1}|U_{2}}(u_{1}|u_{2})=\alpha u_{1}^{2}+(1-\alpha)u_{1}$, $0\leq
u_{1}\leq 1$ (for each $0\leq u_{2}\leq 1$), where
$\alpha\equiv\alpha(u_{2};\theta)=\theta(2u_{2}-1)$. Thus,

$F^{-1}_{U_{1}|U_{2}}(u_{1}|u_{2})=\mathopen{}\mathclose{{}\left\\{\begin{array}[]{c@{\quad\textrm{if}
\ \ }l}u_{1},&\alpha=0,\\\ u_{1,1}\mathds{1}_{\\{0\leq u_{1,1}\leq
1\\}}+u_{1,2}\mathds{1}_{\\{0\leq u_{1,2}\leq 1\\}},&\alpha\neq 0,\\\
\end{array}}\right.$

where $u_{1,1}=[\alpha-1-\\{(1-\alpha)^{2}+4\alpha u_{1}\\}^{1/2}]/(2\alpha)$
and $u_{1,2}=[\alpha-1+\\{(1-\alpha)^{2}+4\alpha u_{1}\\}^{1/2}]/(2\alpha)$.

Let now
$\tilde{\boldsymbol{U}}=(\tilde{U}_{1},\tilde{U}_{2},\tilde{U}_{3})^{T}\sim\tilde{C}_{\theta}$,
see (I.1b). The cdf of $(\tilde{U}_{2},\tilde{U}_{3})^{T}$ is
$F_{\tilde{U}_{2},\tilde{U}_{3}}(u_{2},u_{3})=u_{2}u_{3}$,
$(u_{2},u_{3})^{T}\in[0,1]^{2}$, that is $\tilde{U}_{2}$ and $\tilde{U}_{3}$
are independent. Moreover, the cdf of
$\tilde{U}_{1}|(\tilde{U}_{2},\tilde{U}_{3})^{T}=(u_{2},u_{3})^{T}$ is
$F_{\tilde{U}_{1}|\tilde{U}_{2},\tilde{U}_{3}}(u_{1}|u_{2},u_{3})=\tilde{\alpha}u_{1}^{2}+(1-\tilde{\alpha})u_{1}$,
$0\leq u_{1}\leq 1$ (for each $0\leq u_{2},u_{3}\leq 1$), where
$\tilde{\alpha}\equiv\tilde{\alpha}(u_{2},u_{3};\theta)=-\theta(2u_{2}-1)(2u_{3}-1)$.
Thus,

$\tilde{F}^{-1}_{\tilde{U}_{1}|\tilde{U}_{2},\tilde{U}_{3}}(u_{1}|u_{2},u_{3})=\mathopen{}\mathclose{{}\left\\{\begin{array}[]{c@{\quad\textrm{if}
\ \ }l}u_{1},&\tilde{\alpha}=0,\\\
\tilde{u}_{1,1}\mathds{1}_{\\{0\leq\tilde{u}_{1,1}\leq
1\\}}+\tilde{u}_{1,2}\mathds{1}_{\\{0\leq\tilde{u}_{1,2}\leq
1\\}},&\tilde{\alpha}\neq 0,\\\ \end{array}}\right.$

where
$\tilde{u}_{1,1}=[\tilde{\alpha}-1-\\{(1-\tilde{\alpha})^{2}+4\tilde{\alpha}u_{1}\\}^{1/2}]/(2\tilde{\alpha})$
and
$\tilde{u}_{1,2}=[\tilde{\alpha}-1+\\{(1-\tilde{\alpha})^{2}+4\tilde{\alpha}u_{1}\\}^{1/2}]/(2\tilde{\alpha})$.

From the preceding analysis, if $X_{1}$, $X_{2}$ and $X_{3}$ are independent
$U(0,1)$ distributions, then
$\mathopen{}\mathclose{{}\left(F^{-1}_{U_{1}|U_{2}}(X_{1}|X_{2}),X_{2},X_{3}}\right)^{T}\sim
C_{\theta}$ and
$\mathopen{}\mathclose{{}\left(\tilde{F}^{-1}_{\tilde{U}_{1}|\tilde{U}_{2},\tilde{U}_{3}}(X_{1}|X_{2},X_{3}),X_{2},X_{3}}\right)^{T}\sim\tilde{C}_{\theta}$.

##  Appendix II Bivariate distributions that are used in Subsection 8.2

Bivariate uniform distribution on the circumference of a circle,
$\boldsymbol{U\\{C(0,1)\\}}$

Let $Z_{1}$, $Z_{2}$ be iid standard normal distributions. Consider the rvs
$U_{i}=Z_{i}\big{/}\big{(}Z_{1}^{2}+Z_{2}^{2}\big{)}^{1/2}$, $i=1,2$. Then,
$(U_{1},U_{2})$ is uniformly distributed on the circumference of the circle
$x^{2}+y^{2}=1$ .

Bivariate $\boldsymbol{\exp\\{\lambda_{1},\lambda_{2},\lambda_{12}\\}}$
distribution

Let $\lambda_{1},\lambda_{2},\lambda_{12}>0$ with
$\lambda_{12}<\min\\{\lambda_{1},\lambda_{2}\\}$. Suppose that
$E_{1}\sim\exp(\lambda_{1}-\lambda_{12})$,
$E_{2}\sim\exp(\lambda_{2}-\lambda_{12})$ and $E_{3}\sim\exp(\lambda_{12})$,
and consider the rvs $X_{i}=\min\\{E_{i},E_{3}\\}-1/\lambda_{i}$, $i=1,2$.
Then, $(X_{1},X_{2})$ follows the bivariate exponential distribution.

Bivariate $\textrm{Morgenstern}\boldsymbol{\\{\alpha\\}}$ distribution

Let $\alpha>0$ and $X$, $U$ be independent $U(0,1)$ distributed. Consider the
rvs $Z=\alpha(2X-1)-1$, $W=1-2\alpha(2X-1)+\alpha^{2}(2X-1)^{2}+4\alpha
U(2X-1)$ and $Y=2U/(W^{1/2}-Z)$. Then, $(X,Y)$ follows the bivariate
Morgenstern distribution.

Bivariate $\textrm{Plackett}\boldsymbol{\\{s\\}}$ distribution

Let $s>1$ and $X$, $U$ be independent $U(0,1)$ distributed. Consider the rvs
$W_{1}=U(1-U)$, $W_{2}=s+W_{1}(s-1)^{2}$,
$W_{3}=2W_{1}(s^{2}X+1-X)+s(1-2W_{1})$, $W_{4}=s\\{s+4(1-s)^{2}X(1-X)W_{1}\\}$
and $Y=W_{2}\\{W_{3}-(1-2U)W_{4}^{1/2}\\}/2$. Then, $(X,Y)$ follows the
bivariate Plackett distribution.

Bivariate $\textrm{ALI-HAQ's}\boldsymbol{\\{a,p\\}}$ distribution

Let $a>0$, $0<p<1$ and $X$, $U$ be independent $U(0,1)$ distributed. Consider
the rvs $V=1-U$, $W_{1}=a(2VU+1)+2a^{2}V^{2}U+1$,
$W_{2}=a^{2}(4V^{2}U-4VU+1)+a(4VU-4p+2)+1$ and
$Y=2U(aV-1)^{2}/(W_{1}+W_{2}^{1/2})$. Then, $(X,Y)$ follows the bivariate ALI-
HAQ’s distribution.

Bivariate $\textrm{Gumbel}\boldsymbol{\\{e\\}}$ distribution

Let $e>0$ and $U_{1}\ldots,U_{4}$ be independent $U(0,1)$ distributed.
Consider the rvs $X=-\ln U_{1}$, $W_{1}=1+eX$, $W_{2}=(W_{1}-e)/W_{1}$,
$W_{3}=-\ln U_{2}$ and
$Y=\mathds{1}_{\\{U_{3}<W_{2}\\}}W_{1}W_{3}+\mathds{1}_{\\{U_{3}\geq
W_{2}\\}}W_{1}(W_{3}-\ln U_{4})$. Then, $(X,Y)$ follows the bivariate Gumbel
distribution.

Bivariate $\boldsymbol{t_{5}}$ distribution

Using the package “mvtnorm” in R, we simulate data from bivariate student
distribution with $5$ degree of freedom and variance-covariance matrix
${\mathrel{\mathchoice{\ooalign{$\displaystyle\mathbf{\Sigma}$\cr\hfil$\displaystyle|$\hfil\cr}}{\ooalign{$\textstyle\mathbf{\Sigma}$\cr\hfil$\textstyle|$\hfil\cr}}{\ooalign{$\scriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptstyle|$\hfil\cr}}{\ooalign{$\scriptscriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptscriptstyle|$\hfil\cr}}}}=\big{(}{1\atop
1}~{}{1\atop 4}\big{)}$.

##  Appendix III Numerical results

This sections presents some useful numerical results. Table III.1 contains the
true values of the index $I$ in the normal distribution case
$\mathcal{N}_{3}(\mathbf{0},{\mathrel{\mathchoice{\ooalign{$\displaystyle\mathbf{\Sigma}$\cr\hfil$\displaystyle|$\hfil\cr}}{\ooalign{$\textstyle\mathbf{\Sigma}$\cr\hfil$\textstyle|$\hfil\cr}}{\ooalign{$\scriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptstyle|$\hfil\cr}}{\ooalign{$\scriptscriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptscriptstyle|$\hfil\cr}}}}_{3}(\rho))$,
where the variance/covariance matrix
${\mathrel{\mathchoice{\ooalign{$\displaystyle\mathbf{\Sigma}$\cr\hfil$\displaystyle|$\hfil\cr}}{\ooalign{$\textstyle\mathbf{\Sigma}$\cr\hfil$\textstyle|$\hfil\cr}}{\ooalign{$\scriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptstyle|$\hfil\cr}}{\ooalign{$\scriptscriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptscriptstyle|$\hfil\cr}}}}_{3}(\rho)$
is defined by (2.2b), for various values of $\rho$. Additionally, Figure III.1
presents scatterplots based on random samples from
$\mathcal{N}_{3}(\mathbf{0},{\mathrel{\mathchoice{\ooalign{$\displaystyle\mathbf{\Sigma}$\cr\hfil$\displaystyle|$\hfil\cr}}{\ooalign{$\textstyle\mathbf{\Sigma}$\cr\hfil$\textstyle|$\hfil\cr}}{\ooalign{$\scriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptstyle|$\hfil\cr}}{\ooalign{$\scriptscriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptscriptstyle|$\hfil\cr}}}}_{3}(\rho))$
distributions for various values of $\rho$, giving an illustration of the
shapes of these distributions.

Table III.1: The index $I\equiv I(\rho)$ for various values of $\rho$ when the underline distribution is $\mathcal{N}_{3}\mathopen{}\mathclose{{}\left(\mathbf{0},{\mathrel{\mathchoice{\ooalign{$\displaystyle\mathbf{\Sigma}$\cr\hfil$\displaystyle|$\hfil\cr}}{\ooalign{$\textstyle\mathbf{\Sigma}$\cr\hfil$\textstyle|$\hfil\cr}}{\ooalign{$\scriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptstyle|$\hfil\cr}}{\ooalign{$\scriptscriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptscriptstyle|$\hfil\cr}}}}_{3}(\rho)}\right)$. $\rho\colon$ | 0 | 0.05 | 0.10 | 0.15 | 0.20 | 0.25 | 0.30 | 0.35 | 0.40 | 0.45 | 0.50 | 0.55 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---
$I\colon$ | 0 | 0.027 | 0.052 | 0.077 | 0.101 | 0.124 | 0.148 | 0.171 | 0.194 | 0.218 | 0.243 | 0.267 |
$\rho\colon$ | 0.60 | 0.65 | 0.70 | 0.75 | 0.80 | 0.85 | 0.90 | 0.95 | 0.97 | 0.98 | 0.99 | 0.995 | 1
$I\colon$ | 0.294 | 0.321 | 0.352 | 0.385 | 0.423 | 0.467 | 0.524 | 0.609 | 0.662 | 0.701 | 0.758 | 0.841 | 1

(a) $\rho=0$.

(b) $\rho=0.25$.

(c) $\rho=0.4$.

(d) $\rho=0.5$.

(e) $\rho=0.6$.

(f) $\rho=0.75$.

(g) $\rho=0.8$.

(h) $\rho=1$.

Figure III.1: Scatterplots: Random samples of size $n=1000$ are drawn from the
trivariate normal distributions
$\mathcal{N}_{3}(\mathbf{0},{\mathrel{\mathchoice{\ooalign{$\displaystyle\mathbf{\Sigma}$\cr\hfil$\displaystyle|$\hfil\cr}}{\ooalign{$\textstyle\mathbf{\Sigma}$\cr\hfil$\textstyle|$\hfil\cr}}{\ooalign{$\scriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptstyle|$\hfil\cr}}{\ooalign{$\scriptscriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptscriptstyle|$\hfil\cr}}}}_{3}(\rho))$
for various values of $\rho$.

##  Appendix IV A data set of biomarkers of hepatic injury

Table IV.1: A data set of size $n=208$ of the biomarkers Bilirubin (B), Alanine and Aspartate aminotransferase (ALT and AST) and Alkaline phosphatase (AP). PID | B | AST | ALT | AP |  | PID | B | AST | ALT | AP |  | PID | B | AST | ALT | AP |  | PID | B | AST | ALT | AP
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
1 | $1.1$ | 34 | 54 | 150 |  | 53 | $0.8$ | 29 | 21 | 122 |  | 105 | $0.8$ | 27 | 20 | 63 |  | 157 | $0.4$ | 15 | 16 | 96
2 | $0.8$ | 32 | 23 | 90 |  | 54 | $0.3$ | 26 | 25 | 119 |  | 106 | $0.4$ | 34 | 59 | 143 |  | 158 | $1.1$ | 96 | 75 | 91
3 | $0.4$ | 17 | 21 | 105 |  | 55 | $0.8$ | 59 | 34 | 86 |  | 107 | $4.4$ | 76 | 42 | 126 |  | 159 | $1.2$ | 38 | 28 | 116
4 | $0.2$ | 10 | 17 | 178 |  | 56 | $0.7$ | 171 | 217 | 107 |  | 108 | $3.2$ | 35 | 24 | 118 |  | 160 | $0.5$ | 15 | 18 | 65
5 | $0.5$ | 26 | 38 | 57 |  | 57 | $0.5$ | 15 | 11 | 81 |  | 109 | $0.4$ | 41 | 34 | 47 |  | 161 | $0.2$ | 22 | 19 | 100
6 | $1.6$ | 144 | 116 | 158 |  | 58 | $0.2$ | 13 | 13 | 109 |  | 110 | $1.2$ | 57 | 50 | 115 |  | 162 | $0.3$ | 53 | 47 | 83
7 | $1.1$ | 51 | 28 | 145 |  | 59 | $0.4$ | 51 | 64 | 71 |  | 111 | $1.6$ | 22 | 20 | 63 |  | 163 | $1.6$ | 99 | 56 | 81
8 | $0.5$ | 26 | 15 | 207 |  | 60 | $0.5$ | 15 | 10 | 88 |  | 112 | $2$ | 21 | 18 | 63 |  | 164 | $1.6$ | 37 | 63 | 231
9 | $1.3$ | 48 | 31 | 295 |  | 61 | $0.3$ | 15 | 15 | 121 |  | 113 | $0.5$ | 26 | 38 | 71 |  | 165 | $0.4$ | 40 | 13 | 273
10 | $0.9$ | 25 | 30 | 76 |  | 62 | $1.2$ | 22 | 17 | 61 |  | 114 | $0.2$ | 14 | 15 | 81 |  | 166 | $1$ | 57 | 85 | 189
11 | $0.7$ | 57 | 52 | 88 |  | 63 | $0.8$ | 16 | 18 | 94 |  | 115 | $1.3$ | 35 | 32 | 61 |  | 167 | $6.4$ | 108 | 68 | 108
12 | $0.2$ | 34 | 35 | 76 |  | 64 | $0.2$ | 31 | 30 | 48 |  | 116 | $0.4$ | 15 | 14 | 72 |  | 168 | $1.1$ | 54 | 79 | 91
13 | $0.5$ | 29 | 14 | 66 |  | 65 | $0.5$ | 29 | 25 | 104 |  | 117 | $0.5$ | 30 | 51 | 66 |  | 169 | $0.7$ | 28 | 56 | 202
14 | $1.7$ | 30 | 26 | 121 |  | 66 | $0.6$ | 24 | 18 | 90 |  | 118 | $0.8$ | 40 | 35 | 82 |  | 170 | $0.7$ | 38 | 51 | 60
15 | $1.7$ | 21 | 13 | 94 |  | 67 | $0.3$ | 153 | 183 | 246 |  | 119 | $0.8$ | 46 | 94 | 170 |  | 171 | $0.4$ | 15 | 11 | 106
16 | $1.2$ | 26 | 16 | 207 |  | 68 | $0.8$ | 49 | 121 | 75 |  | 120 | $1$ | 57 | 85 | 189 |  | 172 | $0.4$ | 36 | 45 | 159
17 | $0.9$ | 34 | 38 | 59 |  | 69 | $33.3$ | 121 | 117 | 434 |  | 121 | $1.3$ | 103 | 21 | 71 |  | 173 | $0.6$ | 61 | 47 | 104
18 | $0.8$ | 42 | 89 | 72 |  | 70 | $1$ | 34 | 47 | 68 |  | 122 | $1.5$ | 18 | 22 | 57 |  | 174 | $0.9$ | 18 | 26 | 96
19 | $0.5$ | 50 | 73 | 68 |  | 71 | $0.2$ | 18 | 24 | 48 |  | 123 | $0.6$ | 46 | 61 | 140 |  | 175 | $0.7$ | 49 | 81 | 187
20 | $0.4$ | 43 | 35 | 87 |  | 72 | $0.7$ | 24 | 14 | 97 |  | 124 | $0.7$ | 17 | 19 | 53 |  | 176 | $0.4$ | 22 | 16 | 58
21 | $0.8$ | 19 | 26 | 112 |  | 73 | $4.5$ | 236 | 64 | 242 |  | 125 | $0.8$ | 29 | 30 | 81 |  | 177 | $2.9$ | 51 | 33 | 143
22 | $0.7$ | 16 | 11 | 36 |  | 74 | $0.6$ | 74 | 188 | 72 |  | 126 | $0.7$ | 80 | 93 | 61 |  | 178 | $1.4$ | 24 | 198 | 108
23 | $0.4$ | 20 | 22 | 99 |  | 75 | $0.2$ | 20 | 23 | 80 |  | 127 | $0.7$ | 32 | 22 | 80 |  | 179 | $0.4$ | 29 | 41 | 108
24 | $0.5$ | 40 | 58 | 71 |  | 76 | $0.5$ | 18 | 18 | 58 |  | 128 | $0.8$ | 105 | 170 | 108 |  | 180 | $0.4$ | 26 | 18 | 121
25 | $0.5$ | 44 | 59 | 79 |  | 77 | $0.7$ | 85 | 39 | 89 |  | 129 | $0.6$ | 44 | 50 | 197 |  | 181 | $0.5$ | 24 | 18 | 147
26 | $1.2$ | 30 | 22 | 49 |  | 78 | $1.4$ | 64 | 44 | 183 |  | 130 | $0.3$ | 64 | 24 | 74 |  | 182 | $0.5$ | 81 | 127 | 68
27 | $0.2$ | 20 | 17 | 118 |  | 79 | $1.6$ | 18 | 15 | 73 |  | 131 | $0.6$ | 47 | 45 | 85 |  | 183 | $0.6$ | 26 | 24 | 67
28 | $0.7$ | 36 | 27 | 129 |  | 80 | $0.2$ | 27 | 46 | 175 |  | 132 | $0.4$ | 113 | 77 | 96 |  | 184 | $2.1$ | 55 | 83 | 93
29 | $1.3$ | 58 | 27 | 75 |  | 81 | $1.4$ | 101 | 61 | 44 |  | 133 | $0.9$ | 21 | 16 | 55 |  | 185 | $0.7$ | 26 | 40 | 103
30 | $35.2$ | 182 | 82 | 147 |  | 82 | $0.5$ | 89 | 91 | 125 |  | 134 | $1.1$ | 39 | 26 | 82 |  | 186 | $0.7$ | 23 | 14 | 148
31 | $0.3$ | 25 | 28 | 80 |  | 83 | $1.2$ | 45 | 74 | 91 |  | 135 | $0.7$ | 57 | 69 | 83 |  | 187 | $1.3$ | 21 | 32 | 72
32 | $0.9$ | 20 | 22 | 89 |  | 84 | $1.6$ | 514 | 377 | 876 |  | 136 | $0.9$ | 35 | 56 | 185 |  | 188 | $0.8$ | 34 | 29 | 76
33 | $1.2$ | 55 | 61 | 52 |  | 85 | $0.3$ | 49 | 51 | 104 |  | 137 | $0.7$ | 15 | 30 | 196 |  | 189 | $0.3$ | 56 | 88 | 71
34 | $1.7$ | 34 | 14 | 75 |  | 86 | $0.9$ | 43 | 24 | 70 |  | 138 | $0.3$ | 44 | 49 | 80 |  | 190 | $0.4$ | 23 | 18 | 99
35 | $1.1$ | 35 | 39 | 110 |  | 87 | $0.6$ | 61 | 80 | 59 |  | 139 | $0.6$ | 108 | 133 | 78 |  | 191 | $0.7$ | 39 | 52 | 130
36 | $0.4$ | 21 | 27 | 73 |  | 88 | $0.7$ | 30 | 25 | 68 |  | 140 | $0.6$ | 25 | 19 | 47 |  | 192 | $0.6$ | 51 | 75 | 111
37 | $0.4$ | 24 | 30 | 77 |  | 89 | $0.8$ | 20 | 20 | 71 |  | 141 | $1.2$ | 30 | 23 | 52 |  | 193 | $0.6$ | 21 | 31 | 63
38 | $0.5$ | 38 | 37 | 67 |  | 90 | $0.4$ | 26 | 34 | 57 |  | 142 | $0.6$ | 37 | 41 | 84 |  | 194 | $0.6$ | 25 | 18 | 93
39 | $0.7$ | 22 | 22 | 56 |  | 91 | $0.4$ | 23 | 32 | 147 |  | 143 | $0.3$ | 83 | 76 | 73 |  | 195 | $0.5$ | 82 | 146 | 50
40 | $0.6$ | 53 | 37 | 128 |  | 92 | $0.5$ | 82 | 61 | 67 |  | 144 | $0.6$ | 143 | 274 | 68 |  | 196 | $0.3$ | 67 | 106 | 72
41 | $0.4$ | 25 | 21 | 97 |  | 93 | $0.7$ | 28 | 42 | 64 |  | 145 | $13.4$ | 185 | 82 | 127 |  | 197 | $0.3$ | 19 | 23 | 86
42 | $0.5$ | 52 | 75 | 68 |  | 94 | $2.2$ | 71 | 46 | 181 |  | 146 | $0.4$ | 85 | 183 | 68 |  | 198 | $0.3$ | 24 | 47 | 79
43 | $0.6$ | 28 | 33 | 142 |  | 95 | $0.8$ | 26 | 29 | 95 |  | 147 | $0.6$ | 299 | 190 | 63 |  | 199 | $0.6$ | 46 | 37 | 71
44 | $1$ | 33 | 25 | 59 |  | 96 | $0.5$ | 32 | 39 | 38 |  | 148 | $0.5$ | 35 | 46 | 78 |  | 200 | $0.5$ | 37 | 49 | 94
45 | $0.4$ | 28 | 39 | 164 |  | 97 | $2.2$ | 29 | 33 | 64 |  | 149 | $0.3$ | 14 | 22 | 105 |  | 201 | $0.4$ | 23 | 44 | 69
46 | $0.7$ | 39 | 44 | 84 |  | 98 | $1$ | 40 | 33 | 83 |  | 150 | $0.4$ | 48 | 46 | 69 |  | 202 | $0.5$ | 19 | 20 | 139
47 | $0.8$ | 28 | 41 | 68 |  | 99 | $1.5$ | 52 | 48 | 53 |  | 151 | $0.7$ | 27 | 38 | 99 |  | 203 | $0.9$ | 208 | 346 | 73
48 | $1.2$ | 64 | 26 | 172 |  | 100 | $0.9$ | 85 | 67 | 49 |  | 152 | $0.9$ | 37 | 52 | 50 |  | 204 | $0.6$ | 57 | 64 | 92
49 | $8.4$ | 45 | 47 | 98 |  | 101 | $0.7$ | 30 | 34 | 113 |  | 153 | $0.4$ | 23 | 25 | 70 |  | 205 | $0.5$ | 35 | 60 | 35
50 | $3.7$ | 74 | 102 | 245 |  | 102 | $1.5$ | 16 | 16 | 66 |  | 154 | $0.4$ | 29 | 47 | 99 |  | 206 | $12.7$ | 156 | 134 | 318
51 | $0.8$ | 32 | 35 | 57 |  | 103 | $0.4$ | 51 | 78 | 115 |  | 155 | $3.2$ | 54 | 26 | 110 |  | 207 | $0.7$ | 98 | 139 | 81
52 | $1$ | 293 | 227 | 147 |  | 104 | $0.4$ | 125 | 141 | 140 |  | 156 | $0.7$ | 49 | 88 | 91 |  | 208 | $0.5$ | 15 | 12 | 102

##  Appendix V _R_ codes

The required packages that used.

⬇

rm(list=ls())

library(MASS); require(MASS); sfLibrary(MASS)

library(energy); require(energy)

library(Emcdf); sfLibrary(Emcdf)

library(mvtnorm)

library(gtools)

library(copula)

library("snowfall");library("rlecuyer")

sfInit(parallel=TRUE, cpus=8, type="SOCK")

library(dHSIC)

library(MonoPoly)

The specific formula of $\boldsymbol{1-K_{\varPi}(t)}$ when
$\boldsymbol{d=3}$.

⬇

g <- function(t){ifelse(t==0,1,1-t+t*log(t)-0.5*t*(log(t))^2)}

g <- Vectorize(g)

The pdf $\boldsymbol{f}$ of
$\boldsymbol{\mathcal{N}_{3}(\mathbf{0},{\mathrel{\mathchoice{\ooalign{$\displaystyle\mathbf{\Sigma}$\cr\hfil$\displaystyle|$\hfil\cr}}{\ooalign{$\textstyle\mathbf{\Sigma}$\cr\hfil$\textstyle|$\hfil\cr}}{\ooalign{$\scriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptstyle|$\hfil\cr}}{\ooalign{$\scriptscriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptscriptstyle|$\hfil\cr}}}}_{3}(\boldsymbol{\rho}))}$
distribution.

⬇

f <- function(x,y,z,r12,r13,r23){

S <- matrix(c(1,r12,r13,r12,1,r23,r13,r23,1),3,3)

density <- dmvnorm(c(x,y,z), mean=c(0,0,0), sigma=S)

return(density)}

f <- Vectorize(f)

The cdf $\boldsymbol{F}$ of
$\boldsymbol{\mathcal{N}_{3}(\mathbf{0},{\mathrel{\mathchoice{\ooalign{$\displaystyle\mathbf{\Sigma}$\cr\hfil$\displaystyle|$\hfil\cr}}{\ooalign{$\textstyle\mathbf{\Sigma}$\cr\hfil$\textstyle|$\hfil\cr}}{\ooalign{$\scriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptstyle|$\hfil\cr}}{\ooalign{$\scriptscriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptscriptstyle|$\hfil\cr}}}}_{3}(\boldsymbol{\rho}))}$
distribution.

⬇

F <- function(x,y,z,r12,r13,r23){

S <- matrix(c(1,r12,r13,r12,1,r23,r13,r23,1),3,3)

Pr <- pmvnorm(lower=rep(-Inf, 3), upper=c(x,y,z),

mean=c(0,0,0), corr=S)

return(Pr)}

F <- Vectorize(F)

Code for computing AUK when
$\boldsymbol{\boldsymbol{X}\sim\mathcal{N}_{3}(\mathbf{0},{\mathrel{\mathchoice{\ooalign{$\displaystyle\mathbf{\Sigma}$\cr\hfil$\displaystyle|$\hfil\cr}}{\ooalign{$\textstyle\mathbf{\Sigma}$\cr\hfil$\textstyle|$\hfil\cr}}{\ooalign{$\scriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptstyle|$\hfil\cr}}{\ooalign{$\scriptscriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptscriptstyle|$\hfil\cr}}}}_{3}(\boldsymbol{\rho}))}$.

⬇

AUK <- function(r12,r13,r23){

auk <- integrate(Vectorize(function(x) {

integrate(Vectorize(function(y) {

integrate(function(z)

{g(F(x,y,z,r12,r13,r23))*f(x,y,z,r12,r13,r23)},

-Inf, Inf)$value }), -Inf,Inf)$value }), -Inf,Inf)$value

return(auk)}

Code for computing $\boldsymbol{I}$ when
$\boldsymbol{\boldsymbol{X}\sim\mathcal{N}_{3}\mathopen{}\mathclose{{}\left(\mathbf{0},{\mathrel{\mathchoice{\ooalign{$\displaystyle\mathbf{\Sigma}$\cr\hfil$\displaystyle|$\hfil\cr}}{\ooalign{$\textstyle\mathbf{\Sigma}$\cr\hfil$\textstyle|$\hfil\cr}}{\ooalign{$\scriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptstyle|$\hfil\cr}}{\ooalign{$\scriptscriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptscriptstyle|$\hfil\cr}}}}_{3}(\rho)}\right)}$.

⬇

I <- function(rho){

Sign <- permutations(n=2,r=3,v=c(1,-1), repeats.allowed=TRUE)

auk <- rep(0,8)

for (j in 1:8){

r12 = Sign[j,1]*Sign[j,2]*rho

r13 = Sign[j,1]*Sign[j,3]*rho

r23 = Sign[j,2]*Sign[j,3]*rho

auk[j] <- AUK(r12, r13, r23)}

ind=sqrt(32/57)*sqrt(sum((auk-1/2)^2))

return(ind)}

The ecdf $\boldsymbol{\boldsymbol{\mathbb{F}}_{n}}$ of a data set
$\boldsymbol{\mathbf{X}_{n}=(\boldsymbol{X}_{1},\ldots,\boldsymbol{X}_{n})^{T}}$.

⬇

Fn <- function(Xn){

Fn1 <- function(i){emcdf(Xn, Xn[i,])}

Fn1 <- Vectorize(Fn1)

FFn <- Fn1(c(1:nrow(Xn)))

return(FFn)}

The $\boldsymbol{\hat{\mathrm{AUK}}_{n}}$ estimator based on the data set
$\boldsymbol{\mathbf{X}_{n}=(\boldsymbol{X}_{1},\ldots,\boldsymbol{X}_{n})^{T}}$.

⬇

AUK <- function(Xn){mean(g(Fn(Xn)))}

Code for computing $\boldsymbol{I}$, via Monte Carlo, when
$\boldsymbol{\boldsymbol{X}\sim\mathcal{N}_{3}\mathopen{}\mathclose{{}\left(\mathbf{0},{\mathrel{\mathchoice{\ooalign{$\displaystyle\mathbf{\Sigma}$\cr\hfil$\displaystyle|$\hfil\cr}}{\ooalign{$\textstyle\mathbf{\Sigma}$\cr\hfil$\textstyle|$\hfil\cr}}{\ooalign{$\scriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptstyle|$\hfil\cr}}{\ooalign{$\scriptscriptstyle\mathbf{\Sigma}$\cr\hfil$\scriptscriptstyle|$\hfil\cr}}}}_{3}(\rho)}\right)}$.

⬇

I <- function(N,n,rho){

mu <- c(0,0,0)

Sigma <- rho*matrix(rep(1,9),3,3)+(1-rho)*diag(3)

I1 <- function(i){

if (i==1){

X <- mvrnorm(n, mu, Sigma)

Sign <- permutations(n=2, r=3, v=c(1,-1), repeats.allowed=TRUE)

auk1 <- function(j){AUK(X%*%diag(Sign[j,]))}

auk1 <- Vectorize(auk1)

auk <- auk1(1:8)

A=sqrt(32/57)*sqrt(sum((auk-1/2)^2))

return(A)}}

I1 <- Vectorize(I1)

ind <- mean(I1(rep(1,N)))

return(ind)}

I <- Vectorize(I)

The function $\boldsymbol{\varphi_{3}}$.

⬇

phi <- function(t){1.61*t+4.513*t^2-13.607*t^3+12.235*t^4-3.751*t^5}

phi <- Vectorize(phi)

The $\boldsymbol{\hat{I}_{n}}$ and $\boldsymbol{\hat{I}^{*}_{n}}$ estimators
based on the data set
$\boldsymbol{\mathbf{X}_{n}=(\boldsymbol{X}_{1},\ldots,\boldsymbol{X}_{n})^{T}}$.

⬇

Indices <- function(X){

Sign <- permutations(n=2, r=3, v=c(1,-1), repeats.allowed=TRUE)

auk1 <- function(j){AUK(X%*%diag(Sign[j,]))}

auk1<- Vectorize(auk1)

auk <- auk1(1:8)

ind=sqrt(32/57)*sqrt(sum((auk-1/2)^2))

ind.st <- phi(ind)

return(c(ind,ind.st))}

Estimating the function $\boldsymbol{\varphi_{d}}$ using MonoPoly R package.

⬇

#let varrho and I be as in Step 1 of the algorithm

C <- SOSpol.fit(I, rho, deg.is.odd=TRUE, K=3, a=0, b=1,

monotone = c("increasing", "decreasing"),

trace = FALSE, plot.it = FALSE, type=2,

control = monpol.control())$beta.raw

C <- C-C[1]

C <- C/sum(C)

phi.d <- function(t){evalPol(t,C)}

phi.d <- Vectorize(phi.d)

}

Estimating the standard deviation $\boldsymbol{\sigma_{\varPi}}$.

⬇

Wrapper <- function(d){

l <- length(d)

Fn<- function(Un){

Fn1<- function(i){

emcdf(Un, Un[i,])

}

Fn1<- Vectorize(Fn1)

FFn<- Fn1(c(1:nrow(Un)))

return(FFn)

}

sigma_Pi<- function(r,n,d){

AUK1<- function(i){

if(i==1){

X<- matrix(runif(2*n, min=0, max=1), n, 2)

auk<- mean(pchisq(-2*log(Fn(X)), df=2*d, ncp=0, lower.tail=TRUE))

return(auk)

}

}

AUK1<- Vectorize(AUK1)

AUK<- AUK1(rep(1,r))

Zn<- sqrt(n)*(AUK-1/2)

sigmaPi<- sqrt(var(Zn))

sigma_Pi<- Vectorize(sigma_Pi)

return(sigmaPi)

}

sPi <- sigma_Pi(1000, 50000, d)

}

#

start_time <- Sys.time()

Output<-sfLapply(Value.list, Wrapper)

R <- matrix(unlist(Output), nrow=length(d), ncol=1)

rownames(R) <-d

colnames(R) <-c("sigma_Pi")

end_time <- Sys.time()

running_time<-end_time - start_time

running_time

sfStop()

A code to measure the speed-time of AUK and dHSIC tests

⬇

speed.times <- function(n){

x<- runif(n)

y<- runif(n)

X <- cbind(x,y)

StartTime.AUK <- Sys.time()

AUK <- mean(g(Fn(X)))

Z <- sqrt(n)*abs(AUK-1/2)/(sqrt(19/432))

EndTime.AUK <- Sys.time()

TimeTaken.AUK <- EndTime.AUK - StartTime.AUK

StartTime.dHSIC <- Sys.time()

dHSIC <- dhsic.test(list(x,y), alpha = 0.05, method = "permutation",

kernel = "gaussian", B = 1000, pairwise = FALSE,

bandwidth = 1, matrix.input = FALSE)

EndTime.dHSIC <- Sys.time()

TimeTaken.dHSIC <- EndTime.dHSIC - StartTime.dHSIC

Res=c(n,c(TimeTaken.AUK,TimeTaken.dHSIC))

R <- matrix(Res, nrow=1, ncol=3)

colnames(R) <- c("sample␣size␣n","AUK-time␣(in␣sec)","dHSIC-time␣(in␣sec)")

rownames(R) <- c("")

return(t(R))}