Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Multifunctional Polypyrrole-Based Textile Sensors for Integration into Personal Protection Equipment
Next Article in Special Issue
Artificial Intelligence-Assisted RFID Tag-Integrated Multi-Sensor for Quality Assessment and Sensing
Previous Article in Journal
Progress and Challenge of Sensors for Dairy Food Safety Monitoring
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Distributed MIMO Measurements for Integrated Communication and Sensing in an Industrial Environment

1
Department of Electrical and Information Technology, Lund University, 22100 Lund, Sweden
2
Institute of Communication Networks and Satellite Communications, Graz University of Technology, 8010 Graz, Austria
*
Author to whom correspondence should be addressed.
Sensors 2024, 24(5), 1385; https://doi.org/10.3390/s24051385
Submission received: 31 January 2024 / Revised: 16 February 2024 / Accepted: 19 February 2024 / Published: 21 February 2024
(This article belongs to the Special Issue Sensing Technologies and Wireless Communications for Industrial IoT)

Abstract

:
Many concepts for future generations of wireless communication systems use coherent processing of signals from many distributed antennas. The aim is to improve communication reliability, capacity, and energy efficiency and provide possibilities for new applications through integrated communication and sensing. The large bandwidths available in the higher bands have inspired much work regarding sensing in the millimeter-wave (mmWave) and sub-THz bands; however, the sub-6 GHz cellular bands will still be the main provider of wide cellular coverage due to the more favorable propagation conditions. In this paper, we present a measurement system and results of sub-6 GHz distributed multiple-input-multiple-output (MIMO) measurements performed in an industrial environment. From the measurements, we evaluated the diversity for both large-scale and small-scale fading and characterized the link reliability. We also analyzed the possibility of multistatic sensing and positioning of users in the environment, with the initial results showing a mean-square error below 20  c m on the estimated position. Further, the results clearly showed that new channel models are needed that are spatially consistent and deal with the nonstationary channel properties among the antennas.

1. Introduction

With the advances of the fifth and sixth generation of mobile communication systems, new application fields are emerging such as vehicle-to-everything, machine-to-machine communication, and smart cities [1]. With new applications, requirements on, e.g., data rates, the number of connected devices and reliability increase. Furthermore, many applications need to be able to communicate and sense the environment. In order to do so integrated communication and sensing (ICAS), where the same hardware and spectrum can be used for both purposes, has been given much attention in recent research. An application with very stringent requirements and that could benefit from ICAS is the industrial scenario [2], which is the focus here.
As part of the development of wireless systems, new frequency bands are becoming available for communication, which also enables applications where sensing and communication coexist in the same band and use the same infrastructure and hardware [3,4]. However, due to the more favorable propagation conditions, most systems will probably still operate in the sub-6 GHz band. Furthermore, in fifth-generation networks, massive MIMO technology is seen as the main enabler of many requirements due to its potential to improve signal-to-noise ratio (SNR) and increase coverage due to the array gain, the ability to simultaneously serve many users, and improved reliability; the latter us partly due to the fact that small-scale fading effects effectively vanish as the number of antennas increases. This effect is called channel hardening [5,6].
With the small-scale fading effects significantly reduced, the experienced reliability is to a large extent dependent on the large-scale fading effects. To combat this, distributed (massive) MIMO, where the antennas are spread over a larger physical area, has emerged as a candidate. Solutions such as cell-free massive MIMO [7,8] and holographic MIMO [9] are also being examined. Another candidate is RadioWeaves [3], which is a proposed network architecture that combines distributed arrays and active, large intelligent surfaces [10] with distributed computationto achieve high ultrareliable, low-latency communication. At the same time, the required amount of power to transmit is reduced due to the proximity of the users.
Ultrareliable, low-latency communication is especially important in the industrial scenario. This scenario is a complex and rich environment from a propagation point of view, and channel characterization therefore becomes of great importance for designing radio systems to enable better communication quality and reliability. Hence, fading statistics need to be well studied in a given environment. These fading statistics are of great importance for the design of radio channel models and radio systems and for the development or investigation of, e.g., network schemes and coding techniques [11,12] for a given application.
With large aperture antenna arrays, such as in distributed MIMO, the commonly used assumption of wide-sense stationary and uncorrelated scatterers (WSSUS) channels is no longer valid. There are two types of nonstationarities: (1) The first is related to the large aperture and distributed MIMO, in which the plane wave propagation assumption breaks down and becomes the spherical wave propagation assumption, i.e., an operation in the near field [13]. Different subarrays experience different channels, e.g., due to the various distances to a user and the difference in observing the line-of-sight (LoS) and non-los (NLoS) paths among different antennas. (2) The second is related to the temporal nonstationarity of the environment due to the fact that the channel statistics change over time in dynamic scenarios with the movement of users and other objects. If both of these are violated, the channel is said to be doubly underspread [14].
As new concepts emerge, there is also a need to test the feasibility of these concepts, and an important part of this is designing and building test beds. For distributed MIMO, different designs have been proposed [15,16,17] and channel sounders and/or testbeds have been built. In this work, we take this one step further and present a design of a truly distributed MIMO channel sounder organized as a mesh network where all the links between each antenna in the distributed array can also be measured and exploited for sensing purposes.
With a measurement setup in place, channel measurements are needed to extract the relevant parameters for realistic channel models. In [18,19,20,21,22], measurements of distributed channels were made. For the topic of joint communication and sensing, work has been done mainly in the mmWave bands in [23,24]. Finally, theoretical work and simulations have been performed in [25] for a sub-6 GHz RadioWeave scenario for sensing [26,27]. Most measurements that have been conducted in terms of ICAS have either been performed with a star-shaped design and/or for higher frequencies. Measurements with other topologies and/or sub-6 GHz frequencies are to a large extent lacking.

1.1. Contributions

In this paper, we describe a distributed MIMO channel sounder design. A whole new multi-ink measurement system has been developed to measure the dynamic properties of distributed antenna channels. As in all measurement setups, there is a need for calibration; here, we describe a practical implementation of how this can be done over-the-air (OTA), paving the way for even more advanced OTA calibration algorithms to be developed for more accurate system designs. With this uniquely designed sounder, we conducted a measurement campaign in a realistic and dynamic industrial-like setting. We analyzed the channel characteristics essential for reliability and nonstationarity aspects stemming from the large array and the dynamic environment. Finally, we exploited delay and/or Doppler information in order to explore the possibilities of sub-6   G Hz channels for integrated communication and sensing in a mesh setup.

1.2. Structure of the Paper

In Section 2, the signal model is presented. Then, in Section 3, the developed measurement system is described. We describe the time-division multiple access (TDMA) structure and describe why an automatic gain control (AGC) is implemented. In Section 4, the need for system calibration is discussed as is how the sounder was calibrated in the presented measurement campaign. The measurement campaign is described in Section 5, along with the environment and the channel sounder configuration. The results are presented in Section 6, including an analysis and discussion of both the communication aspect and sensing possibilities of distributed MIMO. Finally, our conclusions and future work are outlined in Section 7.

1.3. Notation

In this paper, [ a ] i and [ A ] i , j denote the ith element of a vector a and the ( i , j ) th entry of a matrix A , respectively. Estimated values are denoted with the hat symbol · ^ . The amplitude of a complex number z is denoted by | z | , z * is the complex conjugate of z, and z is its phase. The Hadamarad product is denoted by ⊙.

2. Signal Model

We consider H a transceiver units distributed in the environment, and their positions are given as p n ( h ) = [ p x , n ( h ) , p y , n ( h ) , p z , n ( h ) ] T R 3 × 1 , with h N a 1 , , H a . In our setup, each transceiver unit supports two independent radio frequency (RF) chains, each connected to a single omnidirectional antenna. It should be noted that a switched—possibly distributed—array can also be connected to the RF chains for even larger setups. In the following signal model, we limit ourselves to the single antenna case for the sake of brevity in notation, but it can easily be extended to the switched array channel sounding system. The H a th unit p n ( H a ) represents the mobile agent, and the other units indexed by h 1 , , H a 1 are the single antenna anchors at known positions. At each time, the h th unit acts as a transmitter and emits a periodic signal s ˜ ( t ) , and the other units p n ( h ) with h N r 1 , , H a h act as receivers. Signals are represented by their complex envelopes with respect to a center frequency f c . The signal received at the h th antenna at the discrete observation time n reads
r n ( h h ) = exp j 2 π μ n ( h h ) t h h exp j 2 π f c ϵ ( h h ) × l = 1 L n α l , n ( h h ) exp j η ( h h ) exp j 2 π f c τ l , n ( h h ) exp j 2 π ν l , n ( h h ) t h h s + w n ( h h ) ,
where the first term comprises L n multipath components (MPCs), l { 1 , , L n } , with each being characterized by its complex amplitude α l , n ( h h ) C and propagation delay τ l , n ( h h ) . Hardware impairments and imperfect synchronization are also characterized in the signal model. More specifically, μ n ( h h ) denotes the frequency offset between the hth unit and the h th unit, η ( h ) denotes the unknown phase offset of the hth unit relative to a reference unit, ϵ ( h ) denotes the time shift due to the clock offset of the hth unit, and ν l , n ( h h ) represents the Doppler shift at the time instant t h h when the channel between the h th transmit antenna and the hth receive antenna of the snapshot n is measured. Note that we omit the frequency dependency of the hardware impairment characteristics, given that a limited signal bandwidth of 40 M Hz is used. Assuming we are transmitting on N f subcarriers, the vector s C N f × 1 accounts for the system response g C N f × 1 and the baseband signal spectrum s f C N f × 1 ; that is, s g s f . The system response is usually measured by a back-to-back calibration procedure. The noise measurement processes w n ( h ) in Equation (1) are independent additive white Gaussian noise (AWGN) with double-sided power spectral density N 0 / 2 .

3. Measurement System

The multilink channel sounder has been developed utilizing the NI-Universal Software Radio Peripheral (USRP) (National Instruments Corporation, Austin, TX, USA) and the software suite LabVIEW 2023. The sounding system is portable and scalable, facilitating various measurement scenarios ranging from indoor and outdoor industrial settings to dense urban environments. The components of our multilink channel sounder system are listed in Table 1 and conceptual overview of the RF parts are shown in Figure 1.
Designed for multilink channel sounding, the channel sounder records and stores all conceivable link combinations between antennas. To avoid interference among links, a TDMA strategy is employed. Each antenna is assigned a unique time slot for signal transmission, during which the remaining antennas are set to receive mode. Figure 2 provides a visual representation of this TDMA structure. As a reference signal, the transmit unit uses a Zadoff–Chu sequence [28]. The signal s ˜ ( t ) is configured as an orthogonal frequency-division multiplexing (OFDM) symbol with Zadoff–Chu samples assigned to subcarriers [15]. The sounding system also allows for the nullification of a specified number of carriers at the spectrum’s edges, thus providing flexibility in bandwidth utilization. The channel sounder captures and streams the raw complex samples directly to the disk on the host computer for subsequent offline postprocessing, which may include symbol averaging.
Tight synchronization is necessary to achieve the TDMA structure. A one-pulse-per-second (1PPS) synchronization signal is distributed to all radios, as well as a stable 10 M Hz reference clock. Depending on the scenario, either synchronized Rubidium clocks or a GPS can be used to discipline the onboard clock and synchronize the triggers. For high-accuracy sensing measurements, rubidium clocks are the preferred choice, we are using SRS FS725 and FS740 (Stanford Research Systems Inc., Sunnyvale, CA, USA).
We record all links, even those that could be assumed to be reciprocal. By saving all channel transfer functions, we enable the possibility to evaluate OTA calibration algorithms.
The sounder is equipped with several adjustable parameters for the TDMA structure, as illustrated in Figure 3. Initially, the reference symbol s ˜ , intended for transmission, is generated and stored in the field-programmable gate array (FPGA) memory of each USRP. Subsequently, the number of repetitions of the reference symbol, denoted with M, is defined. R 2 , with the first symbol effectively serving as a cyclic prefix. Should AGC be used, a description of which will follow, R 3 . This setting accounts for the final symbol’s potential distortion, as hardware adjustments may affect the receive gain during this period. Increasing the value of M can improve the received SNR through symbol averaging. However, this improvement comes at the cost of extended transmission time and a reduced maximum resolvable Doppler frequency.
Furthermore, the structure includes H a TDMA slots, where H a corresponds to the number of antennas (see Figure 2). Following the activation and recording of all elements, the system can enter a quiet state for a duration of B / 120   MHz seconds, where B represents the number of FPGA ticks and 120 M Hz is the FPGA clock rate.

Automatic Gain Control

Figure 4 illustrates three of the H a distributed single antennas. During the first TDMA slot, antenna 1 transmits while all other antennas are in receiving mode. Given the relative distance between antenna 1 and antennas 2 and H a , the latter may require maximization of its receive gain. In the subsequent TDMA slot, the next antenna in the sequence transmits and the rest assume receiving roles. In this TDMA slot, antenna H a , positioned closer to the transmitting antenna, might experience analog to digital converter (ADC) saturation due to the preset gain of the receivers. This elementary example of a realistic scenario illustrates the need for an AGC. Due to the TDMA structure and how the antennas are distributed in space, the gain must be estimated and set within a couple of microseconds. Therefore, the AGC is implemented in the FPGA on board the radio to minimize latency. The implemented AGC is inspired by [17,29].

4. System Calibration

There are several system errors that need to be handled. Some of these errors are more pronounced than are others and stem from different sources, such as temperature variations, clock drift, clock offsets, and timing offset, to mention a few. Many of these errors can be mitigated with well-synchronized measurement equipment and a stable temperature. Now, we briefly describe five different errors and their possible sources. Let us start with the time offset, which simply means that all the USRPs in the system must share the same notion of time so that all saved data from different antennas can be related to each other. If the oscillators on the different transceivers do not provide the exact same carrier frequency, it will result in a carrier frequency offset (CFO). If the aforementioned oscillator phase-locked loops (PLLs) locks in different phases, it will result in a clock phase offset. If the clock frequency of the ADCs is different or imperfect, another frequency offset will be induced, namely, sampling clock frequency offset. Lastly, ADCs might sample the signal at different times due to obtaining the trig signal at different times—e.g., due to different length cables—which is called sampling time offset. A summary can be found in Table 2.
To ensure accurate results, it is essential to perform a back-to-back calibration to remove as much as possible of the described errors. During a back-to-back calibration, the cables from two RF chains are connected together as close as possible to the antenna ports. This gives the transfer function of the complete system between a pair of transceivers. This step must be taken for all combinations of transceiver chains when all radios are operational with the settings intended for use during the measurement campaign. This procedure enables the de-embedding of the radio channel from the measured channel transfer function (CTF), which, in addition to the propagation channel, also includes the influence of cables, connectors, the analog front-end, and digital processing. Furthermore, if possible, it is recommended to characterize the antenna radiation pattern in an anechoic chamber to mitigate the effect of the antenna, thus isolating the wireless CTF. It is important to note that a hardware restart requires a recalibration. This requirement arises due to the reinitialization of the transmit-and-receive chains PLL, which lock onto a random phase after each restart. If the purpose of the sounding is to measure metrics regarding the channel statistics for communication-related evaluation, then this requirement can be relaxed. However, for applications that require the use of coherent signals, such as accurate positioning algorithms, knowing the phase relations between all transceivers is crucial.
The next best thing is to perform a back-to-back calibration after a hardware reset using the same hardware (cables, connectors, etc.) as that used during the measurements. With this method, there will be an error in the phase relationship since it is not possible to control at which phase the PLLs of the different radios will lock. If, for some reason, such as logistical constraints, any back-to-back calibration cannot not be performed, the use of an over-the-air calibration method in postprocessing is required. This approach is feasible if the locations of the antennas are known.
In the following sections, we will describe the steps taken to achieve a calibrated data set. Due to logistical problems, we were unable to perform a back-to-back calibration, and hence we resorted to a combination of post-back-to-back calibration with two of the units to compensate for cable lengths and signal processing time. Then, we applied an over-the-air approach to compensate for the CFO and propagation delay. All steps require that at least a portion of the measured scenario is static so that we can assume that the CTF does not change during the calibration procedure. We ensured that we did not have any time offset by syncing all computers to a local network time protocol (NTP) server. Then, we ensured that the FPGAs shared the same notion of global time. We also assumed that there was no sampling clock frequency offset.

4.1. DC Component

This step is not, in the strict sense, a calibration procedure but is performed because our hardware uses direct down-conversion (DDC) from radio frequency down to baseband. To avoid local oscillator (LO) leakage, the DC component is nulled by transmitting a 0 on the center subchannel ( f c ). Hence, the (complex baseband) DC-component has to be interpolated by taking the average of the amplitude of the two neighboring complex coefficients and the average phase evolution as follows:
r n ( h h ) ^ f = 0 = r n ( h h ) f = Δ f + r n ( h h ) f = Δ f 2 · exp j r n ( h h ) f = Δ f + ϕ ^ n ,
where ( diff a 1 , a 2 , , a L = a 2 a 1 , a 3 a 2 , , a L a L 1 is a function that takes two consecutive values in the vector and their differences, with the resulting vector being one element smaller; unwrap · is the phase unwrapping function of Matlab); and ϕ ^ = E diff unwrap r n ( h h ) is the average phase difference between two consecutive subcarriers after the phase has been unwrapped.

4.2. Carrier Frequency Offset

Even if a good reference clock is provided and distributed, there might be frequency drifts or offsets due, for example, to hardware impairments and/or temperature variations. To identify and remove possible carrier frequency offsets, we use a part of the measurement where all antennas are static. If there is a carrier frequency offset, it comes from the oscillators and is not due to Doppler caused by movements. Inspired by [30], we identify carrier frequency offsets as follows. We collect the snapshots r n ( h h ) received at the antenna h from the antenna h in the N f × N s t matrix H ( h h ) = r 0 ( h h ) , r 1 ( h h ) , r 2 ( h h ) , , r N s t 1 ( h h ) , where n { 0 , 1 , , N st 1 } are all static snapshots, and define the N s t × N s t shift matrix S as follows:
S 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0
This is applied s times to shift the columns of H ( h h ) as follows:
C ( h h ) = H ( h h ) H ( h h ) · S s * .
Discard the s last columns of C ( h h ) since they are all zeros
C s ( h h ) C ( h h ) 1 : N f , 1 : N s t s .
The average carrier frequency offset can now be estimated as follows:
μ ^ ( h h ) = n f = 1 N f n = 1 N s t s C s ( h h ) n f , n
The correction factor then becomes
exp j n μ ^ ( h h ) s , n { 0 , 1 , , N s 1 } .

4.3. Delay Calibration

Assuming line-of-sight conditions with no contributions from MPCs, the transfer function between antennas h and h can be modeled as follows:
r n ( h h ) n f α l , n ( h h ) exp j 2 π f p n ( h ) p n ( h ) c 0 ,
where p n ( h ) p n ( h ) = d n ( h h ) denotes the scalar distance between antennas h and h . Calibrating the delay α ˜ l , n ( h h ) can be omitted.
r n ( h h ) n f = 2 π f · d n ( h h ) c 0 = a · f .
Equation (8) is a straight line with slope a = 2 π d n ( h h ) / c 0 . Since both the frequencies and the constant distance during the snapshots selected for calibration are known, the “true” slope is known. By estimating the measured slope, a ^ , of Equation (8) and with the knowledge of the ground truth positions, a delay calibration can be formulatedas follows:
ϵ ^ ( h h ) a ^ + 2 π d n ( h h ) / c 0 .
Of course, this will not be true in practice, but this is a first approximation to enable calibration to compensate for delays induced by cables and signal processing on the FPGA. If the channel is a non-line-of-sight condition, this approach will overestimate the delay and move the channel impulse response too far.

5. Measurement Campaign

5.1. Environment

The environment for our measurement campaign can be described by a typical industry hall for metal work, e.g., metal lathe, metal cutting. The dimensions are approximately 30 m × 11 m × 8 m (L × W × H), see Figure 5a. There are many metal objects and pieces of machinery that make for a rich wireless environment. Twelve static, frequency coherent, and distributed antennas were divided equally on each long side of the room, approximately 4 m above the floor and separated by 4 m ; see Figure 5b. The infrastructure antennas were tilted approximately 45 degrees to obtain better coverage of the floor level area while strong reflections from the walls directly behind them were avoided. The free-space radiation pattern of the antennas is omnidirectional in cross section, but this will of course not be true as soon as it is attached to the metal structure and other objects close to its proximity. However, as previously mentioned in Section 4, we save and evaluate the radio channel, which is the wireless propagation channel influenced by the antenna radiation pattern. During the measurements, the facility was used as usual, with students working on projects and people moving around.

5.2. Ground Truth

To know where the channel samples are taken and to be able to quantify the accuracy of the radio-based position estimates, a ground truth position is needed. This ground truth usually comes from high-quality global navigation satellite system (GNSS) signals when measurements are performed outdoors. In indoor scenarios, different approaches exist, e.g., camera-based motion capture, use of inertial measurement units (IMU), or a sensor fusion approach using cameras, lasers, and IMUs. However, the acquired ground truth positions must be at least an order of magnitude better than the estimates that are being evaluated.
In our case, all measurements were performed indoors, which ruled out a GNSS solution. Therefore, a combination of a lidar sensor (Ouster OS-Dome 128, Ouster Inc., San Francisco, CA, USA) and an IMU (Microstrain 3DM-GX5-25 (AHRS), Microstrain by HBK, Williston, VT, USA) was used to track the odometry of active and passive users. The sensors were connected to a laptop running Ubuntu 20.04 and the Robot Operating System (ROS) [31] Noetic. All raw sensor messages were saved on disk to allow for the evaluation of different standard SLAM algorithms. In this work, we used the open-source algorithm FAST-LIO2; see [32] for details. The sensors were mounted on a robot to track its position and orientation. During each measured scenario, a new map was constructed as the robot moved. Then, all maps were merged to obtain a single coordinate and reference system. At the mounting point on each of the distributed antennas, we put reflective tape to allow for the easier localization of the antennas on the map; see Figure 6.

5.3. Measured Scenarios

In this study, the moving agent traversed various paths to simulate conditions relevant to robotics or IoT devices in an industrial setting. Each path was traversed multiple times to ensure a robust statistical foundation for channel evaluations. This repetition also facilitates preliminary assessments of data-driven methodologies and machine-learning techniques based on different scenario realizations. All scenarios originated from a fixed position within the environment. For conciseness, this paper focuses on two primary scenarios, henceforth referred to as ref and loop; see Figure 7, which collectively covers critical conditions such as NLoS and LoS links.
In the ref scenario, the robot navigates centrally through the room for approximately 20 s , executes a 180° rotation, and returns to its starting point. The whole sequence lasts around 60 s . In the loop scenario, the robot drives two laps around some machinery and work tables. Parts of the trajectory have a much lower ceiling height than the rest of the hall. The loop scenario lasts approximately 80 s . For both scenarios the parameters of the channel sounder are detailed in Table 3.
We used the remotely controlled robot as an active user in the environment, driving along different paths. The purpose of the measurement was to extract channel statistics for distributed MIMO and to collect channel samples to develop and evaluate positioning algorithms. To achieve this, we collected channel data from routes in the somewhat more open space in the middle of the workshop as well as in the more obstructed parts with blocking from the machinery.

6. Analysis and Discussion

6.1. Maximum Ratio Transmission

To achieve reliable communication with low latency, i.e., no retransmissions, a high SNR is desired. If one has spatial diversity in the form of several transmit and/or receive antennas, as we have in our case, then it is shown that to maximize the SNR at time n, one should use the linear precoder presented in [33]. Collect all uplink snapshots in the matrix
H n = r n ( 1 , 13 ) , , r n ( 12 , 13 ) , C ( N f × H a 1 )
then, using the H a long column vector e = [ 1 , 1 , , 1 ] T consisting of only ones
H n MRT = H n * H n e H n * H n e , C ( N f × 1 )
where the noise is assumed to be white Gaussian and uncorrelated with the signal.
In Figure 8, two representative plots show the channel hardening effect from using distributed antennas; there are no more really deep fading dips. In the ref scenario, we achieved an average array gain of 13.8 dB. In the loop scenario, the average array gain was 14.4 dB. When we averaged all subcarriers, the results were similar, as shown in Figure 9. In the loop scenario, there were still variations in the received power levels in the order of 10 dB, and essentially all antennas experienced NLoS conditions at time 20 s and 58 s . However, despite the challenging propagation conditions, fading levels were small and the received power levels reasonably large. The key takeaway from these results is that there is much to gain from distributing antennas to combat small-scale and large-scale fading which enables reliable communication in challenging environments.

6.2. Local Scattering Function

In the case of a user moving in an industrial environment, the surroundings are usually cluttered with many objects that can have a considerable impact on the behavior of the propagation of wireless signals. Hence, the fading process is nonstationary, which means that the wireless channel can be approximated by a piecewise stationary stochastic process where statistical parameters are valid locally (i.e., in small regions) [14]. To extract parameters from nonstationary channels, we utilize the local scattering function defined in [34,35]. This time-frequency-bounded function covers a stationarity region where the wireless channel can be well approximated by a WSSUS process [14,36].
First, collect all snapshots between the antenna pair ( h h ) in the matrix H ( h h )
H ( h h ) = r 1 ( h h ) , , r N ss ( h h ) , C ( N f × N ss ) ,
Then, following the methodology outlined in [35], we denote the index of the sliding window in time and frequency with k t and k f , respectively. The size of the stationarity region is denoted with M snapshots in the time domain and N samples in the frequency domain. Following [19], we use N = N f and henceforth drop k f . The local scattering function is estimated as in [35]:
C ( h h ) k t ; n , p = 1 I J w = 0 I J 1 H ( h h ) w k t ; n , p 2
where n 0 , , N 1 is the delay index, and p M / 2 , , M / 2 1 is the Doppler index. The local scattering function at k t represents the center value of the time-frequency region. Let H ( h h ) w be the windowed time-variant channel transfer function between each pair of antennas h and h in the stationarity region k t .
H ( h h ) w k t ; n , p = m = M / 2 M / 2 1 q = N / 2 N / 2 1 H ( h h ) m k t , q · G w m , q · e j 2 π p m M n q N ,
where m and q are the relative time and frequency indices within each stationarity region. The relationship between the absolute time index n and the relative time index m is n = k t 1 Δ t + m + M for k t 1 , , N s M Δ t , where Δ t corresponds to the time shift between two consecutive regions of stationarity. The taper functions G w m , q = u i m + M / 2 u ˜ j q + N / 2 are the (separable) band-limited discrete prolate spheroidal sequences (DPSSs) [37], which are well localized within the region M / 2 , M / 2 1 × N / 2 , N / 2 1 . The sequences u i and u ˜ j are indexed by i 0 , , I 1 and j 0 , , J 1 , respectively, with w = i J + j , which is the summation index in Equation (13).
For our measurements, we set I = 1 and J = 2 following the recommendations of [35]. We choose M = 75 , as the region of the minimum-time-stationarity region that corresponds to a duration of 375 m s . Considering the maximum speed of the mobile robot v = 1 m / s , the stationarity region becomes approximately 4.5 wavelengths. As mentioned above, we choose N = N f , assuming that the stationarity bandwidth is 35 M Hz since the relative bandwidth is less than 1%.

6.3. Collinearity

The collinearity metric between the local scattering function in two different time instances allows us investigate the extent of the stationarity region in time, T s n ; that is, how long the WSSUS assumptions will hold [35]. It should be noted that the stationarity time itself will be time dependent due to the changing environment. Stack the N × M elements of C ( h h ) k t ; n , p in a column vector c k t (without the superscript for readability) and define the collinearity R k t 1 , k t 2 as follows:
R k t 1 , k t 2 = c k t 1 T c k t 2 c k t 1 c k t 2 .
As in [35], we define the indicator function γ k t , k ˜ t as
γ k t , k ˜ t = 1 : R k t , k ˜ t > c th , 0 : otherwise .
where a threshold c th is defined. From γ , the (time-varying) stationarity time, T s n , can be estimated as the width of the region around the diagonal. In [35], the threshold c th was somewhat randomly chosen as 0.9. As seen in Figure 10 and Figure 11, we select two (at random) links from the two scenarios ref and loop. We have also plotted how the regions would grow if c th = 0.7 instead.
In scenario ref, the user was moving down in the middle of the workshop, then returning approximately the same path. In Figure 10a,b, we can see that on the way back, we move through a region where the time stationarity region is longer. Here, the channel statistics are valid for a longer distance. In Figure 10c,d, it also looks like the off-diagonal regions indicate that we are actually moving through the same stationarity region on the way back since the collinearity between times 15 s and 45 s is above the threshold.
Performing a similar analysis on the collinearity plots of the loop scenario, where the users performed two complete laps around some machinery, we can also indicate that we are in a stationarity region with similar statistics on the second lap. This becomes more apparent if we lower the threshold, c th , to 0.7; see Figure 11. In general, the stationarity regions seem to become somewhat smaller because of the NLoS conditions.
In Figure 12, we show the corresponding estimated time-varying stationarity regions in meters, T s n , for the two scenarios when c th = 0.7 . The median stationarity region is around 2 m , see Figure 13, which means that the radio channel statistics vary while the UE is moving in the environment. Looking at the recorded statistics and the details of the environment, the rms delay spread, the Doppler power spectrum, and the LoS/NLoS states are changing for just a few meters of movement of the UE, hence the relatively short wide-sense stationarity regions.

6.4. RMS Delay Spread

The power delay profile (PDP), P τ , can be calculated as the marginal expectation of the local scattering function Equation (13) with respect to the Doppler domain as follows:
P ^ τ [ k t ; n ] = 1 M p C ( h h ) k t ; n , p
From this, we can calculate the first and second moments τ and σ τ , respectively, as follows:
σ τ k t = n = 0 N 1 ( n τ s ) 2 P ^ τ [ k t ; n ] n = 0 N 1 P ^ τ [ k t ; n ] τ k t 2 , τ k t = n = 0 N 1 ( n τ s ) P ^ τ [ k t ; n ] n = 0 N 1 P ^ τ [ k t ; n ] ,
where τ s = 1 / ( N Δ f ) . Figure 14 shows the rms delay spread for the different antennas over the two routes, and Figure 15 shows the corresponding CDFs. In calculating the moments in Equation (18), only contributions from the PDP that satisfied certain power thresholds were taken into account [38]. The power threshold was selected as 5 dB above the noise floor to mitigate any spurious component, and 30 dB below the peak to only consider components that had a significant contribution. The median rms delay spread was in the range 38 n s to 54 n s , with significant variations between both antennas and over the routes. We see that the results are in agreement with previous measurements in industry environments [39,40], where the spread was also found to be around 50 n s in a similar-sized environment. In [41], machine-type communication between robot arms was measured in an industry environment. Measurements were made with a bandwidth of 500 M Hz and in a fixed position in the room due to the installation of the robot arm. In the their findings, the delay spread was somewhat lower, around 30 n s . Lastly, in [42], two wideband measurements were performed in what was classified as indoor classroom and industry. The dimensions of the room where the industrial measurements were taken were approximately half those of ours. They reported results of around 70 n s in both LOS and NLOS situations in the industry scenario.

6.5. Doppler Spectral Density

An important metric to characterize dynamic channels is the Doppler spectral density (DSD). In Figure 16, we present the time-variant DSD estimated with MUSIC [43] and ESPRIT [44]. Both methods are so-called super-resolution algorithms and both manage to estimate the Doppler well. There is a model parameter in both algorithms that must be estimate which is related to how many sources (tones) are expected, and this will vary in scenarios such as in the ones presented here. Usually, the model order is estimated using, for example, the Akaike information criterion or minimum description length, but this study, we simply set the model parameter to two. The results showed that even in the challenging parts of the loop scenario, both MUSIC and ESPRIT managed to find the dominant Doppler component.

6.6. Doppler-Delay Positioning and Tracking

To show how well our data set is suited for positioning and to hint at what accuracy can be achieved, we present the initial positioning and tracking results. We focus on an uplink positioning task where the agent transmits signals from which the D-MIMO infrastructure infers its position. For this purpose, we focus on scenario ref, where our aim is to track the agent (×) at a “true” unknown position p n ( 13 ) , moving on a (ground truth) trajectory () based on its uplink signals r n ( h 13 ) received by the D-MIMO antennas () at the “true” known positions { p ( h ) | 1 h 12 } . Before solving the positioning task, we first analyze the data. For this purpose, we collect the available received snapshots until the current step n along the trajectory into overlapping windows with a length of N ν = 150 and assemble them in matrices as follows:
H n ˜ ( h 13 ) = r n N ν + 1 ( h 13 ) , , r n ( h 13 ) C ( N f × N ν ) ,
where we perform a rough delay calibration to account for the time shifts ϵ ( h 13 ) introduced by the clock offsets of the receiving units h w.r.t. the agent. We formulate the ( N f × 1 ) temporal array response in the frequency domain through its elements as follows:
b p n ( 13 ) n f = exp j 2 π f n f τ n ( h 13 ) ,
with f n f denoting n f th frequency bin of the signal in the complex baseband and τ n ( h 13 ) modeling the hypothetical propagation delay from the agent at p n ( 13 ) to the hth receiving antenna at p ( h ) . We further formulate the ( N ν × 1 ) Doppler array response in the time domain through its elements as follows:
c p n ( 13 ) , v n n ˜ = exp j 2 π t n ˜ ν n ( h 13 ) ,
where t n ˜ { 0 , , ( N ν 1 ) / f rep } corresponds to time instances of the current window of snapshots, and ν n ( h 13 ) models the hypothetical Doppler frequency shift depending on the agent position p n ( 13 ) and velocity v n relative to the hth receiving antenna. Note that we omit the dependence on MPCs l > 1 in Equations (20) and (21). In our position and velocity estimator, we model LoS propagation only where NLoS paths enter as disturbance. Since the LoS amplitudes are likely to be large compared to the NLoS amplitudes { α l , n ( h 13 ) | l > 1 } and some of the receiving units h will have the LoS conditions (refer to Figure 5 and Figure 7), the D-MIMO units are likely to jointly agree on the true agent position, even in such an unfavorable industrial environment. We compute the nonphase-coherent empirical Bartlett spectrum (for brevity, we omit the normalization term in the denominator of the classical Bartlett spectrum from Equation (A4)) (see Appendix A for a derivation).
P ^ ( p , v ) = h = 1 12 b H H n ˜ ( h 13 ) c * 2
such that the contributions of all receiving antennas h { 1 , , 12 } are summed noncoherently as powers instead of complex-valued amplitudes because we do not have an accurate phase calibration available between our single-antenna receiving units. In the following, we assume that the agent is moving on a plane parallel to the x y -plane at a known height; hence, we aim for two-dimensional (2D) positioning and velocity estimation in this work, well-knowing that Equation (22) is also suitable for three-dimensional (3D) positioning. We analyze the Bartlett spectrum around an observation step n = 2544 and hence employ a window of received signals n ˜ { 2395 , , 2544 } .
To evaluate the impact of only the temporal array response on the Bartlett spectrum, we choose c : = 1 ( N ν × 1 ) denoting a N ν -dimensional vector of all ones and evaluate Equation (22), which results in the spectrum depicted in Figure 17. Due to the limited bandwidth of B W = 35   M Hz and the respective delay resolution of approximately 8.6   m , the resulting Bartlett spectrum is rather flat around the true agent position (×). Furthermore, imperfections in the delay calibration lead to a bias of the maximum arg max p { P ^ ( p ) } with respect to the true position p n ( 13 ) .
To evaluate the impact of only the Doppler array response on the Bartlett spectrum, we choose b : = 1 ( N f × 1 ) and evaluate Equation (22), resulting in the spectrum shown in Figure 18. At the current time step n, the agent velocity is v n   0.77   m   s 1 . At this speed, the Doppler array response is already much more informative (i.e., it exhibits a higher curvature) around the true agent position than is the temporal array response at the chosen window size. Hence, with a moving agent, the Doppler information quickly dominates over the delay information.
Another benefit of the Doppler array response in Equation (21) is that it also depends on the agent velocity v n , and therefore arg max p , v { P ^ ( p , v ) } is a joint position-velocity estimator. Figure 19 shows the resulting Bartlett spectrum in the velocity domain. At the given speed, the Doppler array response likewise causes a distinct peak around the true agent velocity v n (Δ).
Although reasonable position and velocity estimates can be extracted from a single snapshot of data, state filtering over all snapshots along the trajectory achieves much better results.
To showcase initial results, we employ the empirical Bartlett spectrum from Equation (22) in a particle filter together with a nearly constant velocity state-space model (cf. [45] p. 274). Figure 20 shows the resulting estimatesin the position domain compared to the ground truth trajectory.It is observable that the estimates slowly converge to the true trajectory and follow it closely until the agent performs its 180° turn. In the curve, the velocity decreases, as does the information from the Doppler array response in the position domain because the sensitivity of a Doppler frequency shift with respect to the position; i.e., ν n ( h 13 ) / p n ( 13 ) depends linearly on the (radial) velocity of the agent relative to the hth unit. The estimation accuracy decreases for a moment until the agent moves at maximum speed and the position estimates converge again. Using the Bartlett beamformer-based implementation, we achieve a positioning mean-square error (MSE) of 18.4   c m with respect to our ground truth. These initial results highlight the potential of the dataset for positioning and tracking and set the stage for future work on more elaborate estimators.

7. Conclusions

A new, truly distributed MIMO channel sounding system was developed. The channel sounder was then used to perform measurements in an industry environment. The results show that distributing the antennas will achieve significant channel hardening and avoid deep fading dips due to small-scale and large-scale fading. We also investigated the stationarity regions, in which the WSSUS assumptions held. This showed that the regions are quite small, with stationarity regions in the order of 2 m. We further showed that the RMS delay spread is in line with previous measurements conducted in similar settings and is around 50 n s ; however, it varies between the distributed infrastructure antennas. Also, the Doppler spectral density was investigated by applying two super-resolution algorithms. We showed that in our data, both algorithms can find the dominant Doppler component. Finally, we have highlighted the potential of positioning with D-MIMO in these environments. Despite NLoS conditions, multipath propagation, and rich scattering in an industrial scenario, even a simple Bartlett beamformer can produce promising positioning results with an MSE below 20 c m when paired with a suitable state-space filter. In the future, we will demonstrate a more elaborate estimator that outperforms our current solution. The initial results hint at possible centimeter-level positioning accuracy.
There are several directions for future work from here. Further investigations of channel characteristics are ongoing, where all the link combinations over the measured scenarios are classified as NLoS or LoS and where all available data can contribute to the statistics of the channel. In addition, work investigating the performance of positioning capabilities in the more challenging loop scenario is currently being carried out in parallel to further improvements of the positioning presented in this paper. Another path is to investigate bi- or multistatic radar when the user is device free.

Author Contributions

Conceptualization, C.N. and F.T.; methodology, all; software, C.N., A.F. and B.D.; validation, C.N., A.F. and B.D.; formal analysis, all; investigation, C.N., X.L., A.F. and F.T.; data curation, C.N.; writing—original draft preparation, all; writing—review and editing, all; visualization, C.N. and B.D.; supervision, A.F. and F.T.; project administration, F.T.; funding acquisition, F.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research was partially funded by Connected Systems at Lund University, the strategic research area ELLIIT, and the REINDEER project of the European Union’s Horizon 2020 research and innovation program under grant agreement no. 101013425.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data available on request due to legal restrictions.

Acknowledgments

We would like to thank the Sustainable Production Initiative at Lund University for making their workshop available for the measurements. We also want to thank Associate Anders J. Johansson and Michiel Sandra for their invaluable help during the measurement campaign as well as Sara Willhammar for her discussions, feedback, and proofreading of the material. We also extend our appreciation to Xuesong Cai and the Kungl. Fysiografiska Sällskapet i Lund (Royal Physiographic Society of Lund) for providing funding for the lidar and IMU system.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:
1PPSone pulse per second
AGCautomatic gain control
CPcyclic prefix
D-MIMOdistributed MIMO
DACdigital-to-analog converter
FPGAfield-programmable gate array
ICASintegrated communication and sensing
IMUinertial measurement unit
IoTinternet of things
JCASjoint communication and sensing
lidarlight detection and ranging
LISlarge intelligent surface
LOlocal oscillator
LoSline of sight
LSFlarge-scale fading
M2Mmachine-type communication
mmWavemillimeter wave
MIMOmultiple-input multiple-output
MPCmultipath component
MRCmaximum ratio combining
MRTmaximum ratio transmission
NLoSnon-line of sight
OFDMorthogonal frequency-division multiplexing
OTAover the air
PDFprobability density function
PLpath loss
RbRubidium
RFradio frequency
RMSroot mean square
SLAMsimultaneous localization and mapping
SNRsignal-to-noise ratio
SSFsmall-scale fading
TDMAtime-division multiple access
TOAtime of arrival
UEuser equipment
USRPuniversal software radio peripheral

Appendix A. Doppler-Delay Bartlett Spectrum

Conventionally, the (empirical) Bartlett spectrum is defined as [46]
P ^ ( θ ) = a ( θ ) H R ^ a ( θ ) a H ( θ ) a ( θ ) ,
where a ( θ ) is an array response parameterized on θ (often the angle of arrival of a spatial array response), and
R ^ = 1 N x t = 1 N x x ( t ) x H ( t )
is the sample covariance matrix of N x received signal vectors x .
For multiple parameter estimation, Equation (A1) can be used with a stacked vector of parameterized array responses. In the case of our Doppler-delay Bartlett beamformer, we thus choose a : = vec b ( τ n ) c T ( ν n ) C ( N f N ν × 1 ) and likewise, we stack the received signal matrix into a vector x : = vec H n ˜ C ( N f N ν × 1 ) , where we use only N x = 1 of such vectors to compute Equation (A2). Noting that a H a = a 2 N f N ν , we formulate the Doppler-delay Bartlett spectrum for a single antenna h as follows:
P ^ ( τ n , ν n ) = 1 N f N ν vec b c T H vec H n ˜ vec H n ˜ H vec b c T = 1 N f N ν c b H vec H n ˜ vec H n ˜ H c b = 1 N f N ν b H H n ˜ c * b H H n ˜ c * *
= b H H n ˜ c * 2 N f N ν ,
where we use the identity vec ( A B C ) = C T A vec ( B ) in (A3) with ⊗ denoting the Kronecker product.

References

  1. TSGR. TR 138 913-V15.0.0-5G; Study on Scenarios and Requirements for Next Generation Access Technologies; Technical Report; 3GPP, 2018. [Google Scholar]
  2. Uusitalo, M.A.; Rugeland, P.; Boldi, M.R.; Strinati, E.C.; Demestichas, P.; Ericson, M.; Fettweis, G.P.; Filippou, M.C.; Gati, A.; Hamon, M.H.; et al. 6G Vision, Value, Use Cases and Technologies From European 6G Flagship Project Hexa-X. IEEE Access 2021, 9, 160004–160020. [Google Scholar] [CrossRef]
  3. der Perre, L.V.; Larsson, E.G.; Tufvesson, F.; Strycker, L.D.; Bjornson, E.; Edfors, O. RadioWeaves for efficient connectivity: Analysis and impact of constraints in actual deployments. In Proceedings of the 2019 53rd Asilomar Conference on Signals, Systems, and Computers, Pacific Grove, CA, USA, 3–6 November 2019; pp. 15–22. [Google Scholar] [CrossRef]
  4. Behravan, A.; Baldemair, R.; Parkvall, S.; Dahlman, E.; Yajnanarayana, V.; Bjorkegren, H.; Shrestha, D. Introducing sensing into future wireless communication systems. In Proceedings of the 2022 2nd IEEE International Symposium on Joint Communications & Sensing (JC&S), Seefeld, Austria, 9–10 March 2022; pp. 1–5. [Google Scholar] [CrossRef]
  5. Willhammar, S.; Flordelis, J.; Van der Perre, L.; Tufvesson, F. Channel Hardening in Massive MIMO-A Measurement Based Analysis. In Proceedings of the 2018 IEEE 19th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC), Kalamata, Greece, 25–28 June 2018; pp. 1–5. [Google Scholar] [CrossRef]
  6. Willhammar, S.; Van der Perre, L.; Tufvesson, F. Fading in reflective and heavily shadowed industrial environments with large arrays. arXiv 2022, arXiv:2211.09599. [Google Scholar]
  7. Ngo, H.Q.; Ashikhmin, A.; Yang, H.; Larsson, E.G.; Marzetta, T.L. Cell-Free Massive MIMO Versus Small Cells. IEEE Trans. Wirel. Commun. 2017, 16, 1834–1850. [Google Scholar] [CrossRef]
  8. Björnson, E.; Sanguinetti, L. Scalable Cell-Free Massive MIMO Systems. IEEE Trans. Commun. 2020, 68, 4247–4261. [Google Scholar] [CrossRef]
  9. Huang, C.; Hu, S.; Alexandropoulos, G.C.; Zappone, A.; Yuen, C.; Zhang, R.; Renzo, M.D.; Debbah, M. Holographic MIMO Surfaces for 6G Wireless Networks: Opportunities, Challenges, and Trends. IEEE Wirel. Commun. 2020, 27, 118–125. [Google Scholar] [CrossRef]
  10. Hu, S.; Rusek, F.; Edfors, O. Beyond Massive MIMO: The Potential of Data Transmission with Large Intelligent Surfaces. IEEE Trans. Signal Process. 2018, 66, 2746–2758. [Google Scholar] [CrossRef]
  11. Oestges, C. Multi-link propagation modeling for beyond next generation wireless. In Proceedings of the 2011 Loughborough Antennas & Propagation Conference, Loughborough, UK, 14–15 November 2011; pp. 1–8. [Google Scholar]
  12. Wilding, T.; Deutschmann, B.J.B.; Nelson, C.; Li, X.; Tufvesson, F.; Witrisal, K. Propagation Modeling for Physically Large Arrays: Measurements and Multipath Component Visibility. In Proceedings of the 2023 Joint European Conference on Networks and Communications & 6G Summit (EuCNC/6G Summit), Gothenburg, Sweden, 6–9 June 2023; pp. 1–6. [Google Scholar]
  13. Fedorov, A.; Zhang, H.; Chen, Y. User Localization Using Random Access Channel Signals in LTE Networks with Massive MIMO. In Proceedings of the 2018 27th International Conference on Computer Communication and Networks (ICCCN), Hangzhou, China, 30 July–2 August 2018; pp. 1–9. [Google Scholar] [CrossRef]
  14. Matz, G. Doubly underspread non-WSSUS channels: Analysis and estimation of channel statistics. In Proceedings of the 2003 4th IEEE Workshop on Signal Processing Advances in Wireless Communications-SPAWC 2003 (IEEE Cat. No.03EX689), Rome, Italy, 15–18 June 2003. [Google Scholar] [CrossRef]
  15. Wassie, D.A.; Rodriguez, I.; Berardinelli, G.; Tavares, F.M.; Sørensen, T.B.; Hansen, T.L.; Mogensen, P. An agile multi-node multi-antenna wireless channel sounding system. IEEE Access 2019, 7, 17503–17516. [Google Scholar] [CrossRef]
  16. Zelenbaba, S.; Löschenbrand, D.; Hofer, M.; Dakić, A.; Rainer, B.; Humer, G.; Zemen, T. A scalable mobile multi-node channel sounder. In Proceedings of the IEEE Wireless Commun. and Networking Conference (WCNC) (WCNC), Seoul, Republic of Korea, 25–28 May 2020; pp. 1–6. [Google Scholar] [CrossRef]
  17. Stanko, D.; Sommerkorn, G.; Ihlow, A.; Del Galdo, G. Enable Software-Defined Radios for Real-Time MIMO Channel Sounding. In Proceedings of the 2021 IEEE International Instrumentation and Measurement Technology Conference (I2MTC), Glasgow, UK, 17–20 May 2021; pp. 1–5. [Google Scholar] [CrossRef]
  18. Löschenbrand, D.; Hofer, M.; Bernadó, L.; Humer, G.; Schrenk, B.; Zelenbaba, S.; Zemen, T. Distributed Massive MIMO Channel Measurements in Urban Vehicular Scenario. In Proceedings of the 2019 13th European Conference on Antennas and Propagation (EuCAP), Krakow, Poland, 31 March–5 April 2019; pp. 1–5. [Google Scholar]
  19. Zelenbaba, S.; Rainer, B.; Hofer, M.; Löschenbrand, D.; Dakić, A.; Bernadó, L.; Zemen, T. Multi-Node Vehicular Wireless Channels: Measurements, Large Vehicle Modeling, and Hardware-in-the-Loop Evaluation. IEEE Access 2021, 9, 112439–112453. [Google Scholar] [CrossRef]
  20. Guevara, A.P.; Bast, S.D.; Pollin, S. Weave and Conquer: A Measurement-based Analysis of Dense Antenna Deployments. In Proceedings of the ICC 2021—IEEE International Conference on Communications, Montreal, QC, Canada, 14–23 June 2021. [Google Scholar] [CrossRef]
  21. Löschenbrand, D.; Hofer, M.; Bernadó, L.; Zelenbaba, S.; Zemen, T. Towards Cell-Free Massive MIMO: A Measurement-Based Analysis. IEEE Access 2022, 10, 89232–89247. [Google Scholar] [CrossRef]
  22. Nelson, C.; Li, X.; Wilding, T.; Deutschmann, B.; Witrisal, K.; Tufvesson, F. Large Intelligent Surface Measurements for Joint Communication and Sensing. In Proceedings of the 2023 Joint European Conference on Networks and Communications & 6G Summit (EuCNC/6G Summit), Gothenburg, Sweden, 6–9 June 2023; pp. 228–233. [Google Scholar] [CrossRef]
  23. Sankar, R.S.P.; Deepak, B.; Chepuri, S.P. Joint Communication and Radar Sensing with Reconfigurable Intelligent Surfaces. In Proceedings of the IEEE 22nd Int. Workshop on Sig. Proc. Adv. in Wireless Comm. (SPAWC), Lucca, Italy, 27–30 September 2021; pp. 471–475. [Google Scholar] [CrossRef]
  24. Zhang, J.A.; Rahman, M.L.; Wu, K.; Huang, X.; Guo, Y.J.; Chen, S.; Yuan, J. Enabling Joint Communication and Radar Sensing in Mobile Networks—A Survey. IEEE Commun. Surv. Tutor. 2022, 24, 306–345. [Google Scholar] [CrossRef]
  25. Fascista, A.; Deutschmann, B.J.; Keskin, M.F.; Wilding, T.; Coluccia, A.; Witrisal, K.; Leitinger, E.; Seco-Granados, G.; Wymeersch, H. Uplink Joint Positioning and Synchronization in Cell-Free Deployments with Radio Stripes. arXiv 2023, arXiv:2302.03387. [Google Scholar]
  26. Fouda, A.; Keating, R.; Cha, H.S. Toward cm-Level Accuracy: Carrier Phase Positioning for IIoT in 5G-Advanced NR Networks. In Proceedings of the 2022 IEEE 33rd Annual International Symposium on Personal, Indoor and Mobile Radio Communications (PIMRC), Kyoto, Japan, 12–15 September 2022; pp. 782–787. [Google Scholar] [CrossRef]
  27. Wymeersch, H.; Amiri, R.; Seco-Granados, G. Fundamental Performance Bounds for Carrier Phase Positioning in Cellular Networks. arXiv 2023, arXiv:2306.12133. [Google Scholar]
  28. Chu, D. Polyphase codes with good periodic correlation properties (Corresp.). IEEE Trans. Inf. Theory 1972, 18, 531–532. [Google Scholar] [CrossRef]
  29. Sobaihi, K.; Hammoudeh, A.; Scammell, D. Automatic gain control on FPGA for software-defined radios. In Proceedings of the Wireless Telecommunications Symposium 2012, London, UK, 18–20 April 2012; pp. 1–4. [Google Scholar] [CrossRef]
  30. Tufvesson, F.; Edfors, O.; Faulkner, M. Time and frequency synchronization for OFDM using PN-sequence preambles. In Proceedings of the Gateway to 21st Century Communications Village. VTC 1999-Fall. IEEE VTS 50th Vehicular Technology Conference (Cat. No.99CH36324), Amsterdam, The Netherlands, 19–22 September 1999; Volume 4, pp. 2203–2207. [Google Scholar] [CrossRef]
  31. Quigley, M.; Conley, K.; Gerkey, B.; Faust, J.; Foote, T.; Leibs, J.; Wheeler, R.; Ng, A.Y. ROS: An open-source Robot Operating System. In Proceedings of the ICRA Workshop on Open Source Software, Kobe, Japan, 12–17 May 2009; Volume 3, p. 5. [Google Scholar]
  32. Xu, W.; Cai, Y.; He, D.; Lin, J.; Zhang, F. FAST-LIO2: Fast Direct LiDAR-Inertial Odometry. IEEE Trans. Robot. 2022, 38, 2053–2073. [Google Scholar] [CrossRef]
  33. Lo, T. Maximum ratio transmission. IEEE Trans. Commun. 1999, 47, 1458–1461. [Google Scholar] [CrossRef]
  34. Matz, G. On non-WSSUS wireless fading channels. IEEE Trans. Wirel. Commun. 2005, 4, 2465–2478. [Google Scholar] [CrossRef]
  35. Bernadó, L. Non-Stationarity in Vehicular Wireless Channels. Ph.D. Thesis, Technische Universität Wien, Vienna, Austria, 2012. [Google Scholar]
  36. Matz, G. Characterization of non-WSSUS fading dispersive channels. In Proceedings of the IEEE International Conference on Communications, 2003 ICC’03, Anchorage, AK, USA, 11–15 May 2003. [Google Scholar] [CrossRef]
  37. Slepian, D. Prolate spheroidal wave functions, fourier analysis, and uncertainty—V: The discrete case. Bell Syst. Tech. J. 1978, 57, 1371–1430. [Google Scholar] [CrossRef]
  38. Czink, N. The Random-Cluster Model: A Stochastic MIMO Channel Model for Broadband Wireless Communication Systems of the 3rd Generation and Beyond. Ph.D. Thesis, Technische Universität Wien, Vienna, Austria, 2007. [Google Scholar]
  39. Rappaport, T.S.; Seidel, S.Y.; Takamizawa, K. Statistical channel impulse response models for factory and open plan building radio communicate system design. IEEE Trans. Commun. 1991, 39, 794–807. [Google Scholar] [CrossRef]
  40. Karedal, J.; Wyne, S.; Almers, P.; Tufvesson, F.; Molisch, A.F. UWB channel measurements in an industrial environment. In Proceedings of the IEEE Global Telecommunications Conference, 2004 GLOBECOM’04, Dallas, TX, USA, 29 November–3 December 2004; Volume 6, pp. 3511–3516. [Google Scholar]
  41. Holfeld, B.; Wieruch, D.; Raschkowski, L.; Wirth, T.; Pallasch, C.; Herfs, W.; Brecher, C. Radio channel characterization at 5.85 GHz for wireless M2M communication of industrial robots. In Proceedings of the 2016 IEEE Wireless Communications and Networking Conference, Doha, Qatar, 3–6 April 2016; pp. 1–7. [Google Scholar] [CrossRef]
  42. Zhang, Q.; Loh, T.H.; Zhou, D.; Shen, F.; Huang, Z.; Qin, F. A Measurement Campaign of Industrial Environment for Ultra Reliable IIoT Systems. In Proceedings of the 2022 14th International Conference on Wireless Communications and Signal Processing (WCSP), Nanjing, China, 1–3 November 2022; pp. 644–648. [Google Scholar] [CrossRef]
  43. Schmidt, R. Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas Propag. 1986, 34, 276–280. [Google Scholar] [CrossRef]
  44. Paulraj, A.; Roy, R.; Kailath, T. Estimation Of Signal Parameters Via Rotational Invariance Techniques- Esprit. In Proceedings of the Nineteeth Asilomar Conference on Circuits, Systems and Computers, Pacific Grove, CA, USA, 6–8 November 1985; pp. 83–89. [Google Scholar] [CrossRef]
  45. Bar-Shalom, Y.; Li, X.-R.; Kirubarajan, T. Estimation with Applications to Tracking and Navigation; Wiley: Hoboken, NJ, USA, 2002. [Google Scholar] [CrossRef]
  46. Krim, H.; Viberg, M. Two decades of array signal processing research: The parametric approach. IEEE Signal Process. Mag. 1996, 13, 67–94. [Google Scholar] [CrossRef]
Figure 1. An illustration of a three-node, multilink setup. The dashed line between the two Rubidium clocks (Rb-clock 1 and Rb-clock 2) illustrates that if the two Rubidium clocks are well synchronized—over several hours—then they can be disconnected for some time without losing the synchronization of the radios. To the RF ports of the USRPs, one can either connect single antennas or switched arrays.
Figure 1. An illustration of a three-node, multilink setup. The dashed line between the two Rubidium clocks (Rb-clock 1 and Rb-clock 2) illustrates that if the two Rubidium clocks are well synchronized—over several hours—then they can be disconnected for some time without losing the synchronization of the radios. To the RF ports of the USRPs, one can either connect single antennas or switched arrays.
Sensors 24 01385 g001
Figure 2. During one TDMA slot, only one antenna is transmitting while all others are receiving. In the next TDMA slot, the next antenna is transmitting while all other are listening. By saving all channels, even the reciprocal ones, one can use the information for over-the-air calibration.
Figure 2. During one TDMA slot, only one antenna is transmitting while all others are receiving. In the next TDMA slot, the next antenna is transmitting while all other are listening. By saving all channels, even the reciprocal ones, one can use the information for over-the-air calibration.
Sensors 24 01385 g002
Figure 3. The TDMA-based signal structure. Each antenna is assigned a dedicated TDMA slot. During each transmission, the antenna transmits R repetitions of the sounding signal, with some being used as guards and the rest for averaging to increase the signal-to-noise ratio.
Figure 3. The TDMA-based signal structure. Each antenna is assigned a dedicated TDMA slot. During each transmission, the antenna transmits R repetitions of the sounding signal, with some being used as guards and the rest for averaging to increase the signal-to-noise ratio.
Sensors 24 01385 g003
Figure 4. Depiction of three of the H a antennas distributed in space. In TDMA slot 1, the antenna 1 is transmitting while all the others are listening. In then next TDMA slot, antenna 2 is transmitting. Since antennas are distributed in space, it is clear from the figure that an AGC is needed; antenna H a might need all the available gain when antenna 1 is transmitting, while that same gain setting might saturate the ADC when antenna 2 is transmitting.
Figure 4. Depiction of three of the H a antennas distributed in space. In TDMA slot 1, the antenna 1 is transmitting while all the others are listening. In then next TDMA slot, antenna 2 is transmitting. Since antennas are distributed in space, it is clear from the figure that an AGC is needed; antenna H a might need all the available gain when antenna 1 is transmitting, while that same gain setting might saturate the ADC when antenna 2 is transmitting.
Sensors 24 01385 g004
Figure 5. (a) A photograph showing a view of the environment. The hall is approximately 30 m × 11  m with a ceiling height between 8  m   to 10  m depending on location. (b) A photograph showing the placement of four antennas (circled in red). In total, there were twelve distributed antennas; six on each side of the hall. They were situated 4  m above the floor, with a separation of 4  m .
Figure 5. (a) A photograph showing a view of the environment. The hall is approximately 30 m × 11  m with a ceiling height between 8  m   to 10  m depending on location. (b) A photograph showing the placement of four antennas (circled in red). In total, there were twelve distributed antennas; six on each side of the hall. They were situated 4  m above the floor, with a separation of 4  m .
Sensors 24 01385 g005
Figure 6. Example output from FAST-LIO2. The figure depicts the intensity of the points in the scan, and the extra intensity spots on antenna locations with the help of reflective tape, as seen in the circled regions. The cyan colored path shown is an example output of the ground truth position of the user.
Figure 6. Example output from FAST-LIO2. The figure depicts the intensity of the points in the scan, and the extra intensity spots on antenna locations with the help of reflective tape, as seen in the circled regions. The cyan colored path shown is an example output of the ground truth position of the user.
Sensors 24 01385 g006
Figure 7. A top-down schematic of the workshop where the measurements were performed. The 12 static antennas are labeled according to the figure. The two paths are depicted where the samples were taken. Machinery or equipment are marked purple, except the machine colored in green which is tall enough to put antenna 8 in NLoS for the majority of the measurements. The part colored in gray is where the ceiling is much lower than the rest of the workshop.
Figure 7. A top-down schematic of the workshop where the measurements were performed. The 12 static antennas are labeled according to the figure. The two paths are depicted where the samples were taken. Machinery or equipment are marked purple, except the machine colored in green which is tall enough to put antenna 8 in NLoS for the majority of the measurements. The part colored in gray is where the ceiling is much lower than the rest of the workshop.
Sensors 24 01385 g007
Figure 8. Using maximum ratio transmission (MRT), we clearly see a channel hardening effect where the deep fading dips are canceled. For visualization, a subset of the snapshots from the ref scenario is selected, and only one, antenna 8 is plotted for comparison with MRT. To the left is MRT for a section of scenario ref, and on the right, scenario loop is shown.
Figure 8. Using maximum ratio transmission (MRT), we clearly see a channel hardening effect where the deep fading dips are canceled. For visualization, a subset of the snapshots from the ref scenario is selected, and only one, antenna 8 is plotted for comparison with MRT. To the left is MRT for a section of scenario ref, and on the right, scenario loop is shown.
Sensors 24 01385 g008
Figure 9. Using MRT, we clearly see a channel hardening effect where the deep fading dips are canceled. Here, all links are depicted using the average power of all subcarriers.
Figure 9. Using MRT, we clearly see a channel hardening effect where the deep fading dips are canceled. Here, all links are depicted using the average power of all subcarriers.
Sensors 24 01385 g009
Figure 10. The collinearity for two different links of the ref scenario, with different threshold values. (a) Link (3, 13) with threshold 0.9, (b) link (3, 13) with threshold 0.7, (c) link (9, 13) with threshold 0.7, and (d) link (9, 13) with threshold 0.7.
Figure 10. The collinearity for two different links of the ref scenario, with different threshold values. (a) Link (3, 13) with threshold 0.9, (b) link (3, 13) with threshold 0.7, (c) link (9, 13) with threshold 0.7, and (d) link (9, 13) with threshold 0.7.
Sensors 24 01385 g010aSensors 24 01385 g010b
Figure 11. The collinearity for two different links of the loop scenario, with different threshold values. (a) Link (3, 13) with threshold 0.9, (b) link (3, 13) with threshold 0.7, (c) link (9, 13) with threshold 0.7, and (d) link (9, 13) with threshold 0.7.
Figure 11. The collinearity for two different links of the loop scenario, with different threshold values. (a) Link (3, 13) with threshold 0.9, (b) link (3, 13) with threshold 0.7, (c) link (9, 13) with threshold 0.7, and (d) link (9, 13) with threshold 0.7.
Sensors 24 01385 g011
Figure 12. The time-varying stationarity region (in m) with a threshold of 0.7.
Figure 12. The time-varying stationarity region (in m) with a threshold of 0.7.
Sensors 24 01385 g012
Figure 13. The CDF of the time-varying stationarity region (in m) with a threshold of 0.7.
Figure 13. The CDF of the time-varying stationarity region (in m) with a threshold of 0.7.
Sensors 24 01385 g013
Figure 14. RMS Delay spread for the two scenarios, calculated using the local scattering function.
Figure 14. RMS Delay spread for the two scenarios, calculated using the local scattering function.
Sensors 24 01385 g014
Figure 15. The empirical CDF for the RMS delay spread for the two scenarios, calculated using the local scattering function.
Figure 15. The empirical CDF for the RMS delay spread for the two scenarios, calculated using the local scattering function.
Sensors 24 01385 g015
Figure 16. The normalized Doppler spectral density (DSD) estimated with MUSIC for link (4, 13). The red line is the theoretical LoS Doppler.
Figure 16. The normalized Doppler spectral density (DSD) estimated with MUSIC for link (4, 13). The red line is the theoretical LoS Doppler.
Sensors 24 01385 g016
Figure 17. Bartlett spectrum in the position domain exploiting only delay information.
Figure 17. Bartlett spectrum in the position domain exploiting only delay information.
Sensors 24 01385 g017
Figure 18. Bartlett spectrum in the position domain exploiting only Doppler information.
Figure 18. Bartlett spectrum in the position domain exploiting only Doppler information.
Sensors 24 01385 g018
Figure 19. Bartlett spectrum in the velocity domain exploiting only Doppler information.
Figure 19. Bartlett spectrum in the velocity domain exploiting only Doppler information.
Sensors 24 01385 g019
Figure 20. Initial trajectory estimation result based on the Bartlett estimator using Equation (22) with a particle filter. The particle filter is initialized at a position p = [ 9.5 , 6.5 ] T , close to its first estimate p ^ 1 ( 13 ) .
Figure 20. Initial trajectory estimation result based on the Bartlett estimator using Equation (22) with a particle filter. The particle filter is initialized at a position p = [ 9.5 , 6.5 ] T , close to its first estimate p ^ 1 ( 13 ) .
Sensors 24 01385 g020
Table 1. Hardware for the multilink measurement system.
Table 1. Hardware for the multilink measurement system.
HardwareAmountDescription
NI-USRP 2953r 40  M Hz  (National Instruments Corporation, Austin, TX, USA)7USRP
SRS FS725 (Stanford Research Systems Inc., Sunnyvale, CA, USA)310  M Hz and 1 pulse per second (1PPS) Rb standard
SRS FS740 (Stanford Research Systems Inc., Sunnyvale, CA, USA)110  M Hz and 1PPS with GNSS
Host computers7Radio control and logging data
Hoverboard1Acting as mobile agent/mobile user (UE)
Joymax SAF-6571RS3X antennas (Joymax Electronics Co., Ltd., Tao-yuan City, Taiwan)1312 as infrastructure and 1 on the UE
Ouster OSDome (128 lines) (Ouster Inc., San Francisco, CA, USA)1The light detection and ranging (lidar) used for simultaneous localization and mapping (SLAM)
Microstrain 3DM-GX5-25 (AHRS) (Microstrain by HBK, Williston, VT, USA)19-DoF IMU for SLAM
Table 2. Summary of measurement errors.
Table 2. Summary of measurement errors.
System ErrorsSource
Carrier Frequency Offset (CFO)The oscillators do not provide the same frequency.
Clock Phase OffsetThe PLLs lock on random—and different—phases.
Sampling Clock frequency OffsetThe clock frequency of the ADCs are not the same.
Sampling Time OffsetThe ADCs samples at different times.
Time OffsetThe system does not share the same notion of time.
Table 3. Channel sounding parameters.
Table 3. Channel sounding parameters.
Parameter DescriptionValueParameter DescriptionValue
Number of antennas, H a 13Carrier frequency, f c 3.75   G Hz
Frequency spacing, Δ f 78.125   k Hz Bandwidth, B W 40  M Hz
Active subcarriers, N f 449Number of subcarriers, N sc 512
Signal length, τ max 12.8   μ s Signal repetitions, R4
Snapshot length, H a · R · τ max 665.6   μ s Repetition rate, f rep 200  Hz  (5  m s )
Max. resolvable velocity, v max m / s Transmit power, P TX 19  dBm
Measurement length, T T 60 , 80 Signal spacing, quite 4334.4   μ s
Digital-to-analog back-off, A DAC 0.9
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Nelson, C.; Li, X.; Fedorov, A.; Deutschmann, B.; Tufvesson, F. Distributed MIMO Measurements for Integrated Communication and Sensing in an Industrial Environment. Sensors 2024, 24, 1385. https://doi.org/10.3390/s24051385

AMA Style

Nelson C, Li X, Fedorov A, Deutschmann B, Tufvesson F. Distributed MIMO Measurements for Integrated Communication and Sensing in an Industrial Environment. Sensors. 2024; 24(5):1385. https://doi.org/10.3390/s24051385

Chicago/Turabian Style

Nelson, Christian, Xuhong Li, Aleksei Fedorov, Benjamin Deutschmann, and Fredrik Tufvesson. 2024. "Distributed MIMO Measurements for Integrated Communication and Sensing in an Industrial Environment" Sensors 24, no. 5: 1385. https://doi.org/10.3390/s24051385

APA Style

Nelson, C., Li, X., Fedorov, A., Deutschmann, B., & Tufvesson, F. (2024). Distributed MIMO Measurements for Integrated Communication and Sensing in an Industrial Environment. Sensors, 24(5), 1385. https://doi.org/10.3390/s24051385

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop