Measurement Based Stochastic Channel Model for 60 GHz Mmwave Industrial Communications

Communications in the mmWave spectrum are gaining relevance in the last years as they are a promising candidate to cope with the increasing demand of throughput and latency in different use cases. Nowadays, several efforts have been carried out to characterize the propagation medium of these signals with the aim of designing their corresponding communication protocols accordingly, and a wide variety of both outdoor/indoor locations have already been studied. However, very few works endorse industrial scenarios, which are particularly demanding due to their stringent requirements in terms of reliability, determinism, and latency. This work aims to provide an insight of the propagation of 60 GHz mmWave signals in a typical industrial workshop in order to explore the particularities of this kind of scenario. In order to achieve this, an extensive measurement campaign has been carried out in this environment and a stochastic channel model has been proposed and validated.


I. INTRODUCTION
Communications in the mmWave spectrum  are gaining notoriety in the last years, as they are a promising candidate to cope with the increasing demand in different sectors as mobile internet, multimedia, 5G, autonomous driving, smart manufacturing, virtual reality, and healthcare, among others [1].The main advantage of these frequencies is the huge amount of available license-free bandwidth, which allows to increase the throughput and robustness of communications.Several standards that operate over these frequencies can be found in the literature [2], being the most extended ones 5G New Radio (NR) [3] and 802.11 (amendments 802.11ad [4] and 802.11ay [5]).
An important aspect for the development of standards as the ones mentioned above is the propagation medium, also known as the channel.A characterization and modeling of this aspect is necessary to estimate the intensity of small and large scale fading effects (shadowing, multipath, etc.) in order to design and validate the protocols accordingly.This process can be carried out in two different manners.The first alternative uses the propagation properties of materials and the fundamental laws of electromagnetism for the modeling process leading to deterministic results.Techniques as ray tracing fall within this category.The other method uses on site measurements to extract metrics in any given scenario, which are later used to adjust parameters of stochastic channel models.Some examples of this kind would be the 3GPP [6] and WINNER II [7] models.
Nowadays, several measurement campaigns and modeling efforts can be found in the literature for the mmWave spectrum.These works cover a wide variety of scenarios for the different industrial, scientific and medical (ISM) bands, where both outdoor [8] (street, rural, campus, etc.) and indoor [9] (residence, office, etc.) sites are considered.Some works present measurements in industrial scenarios but these are rather scarce, as will be described later in this section.
The different use cases for communication systems in industry are oriented to managing the processes that are carried out within a factory [10].This includes tasks as monitoring and remote control of elements as production lines, robotic arms or logistic robots, and the communications related to them need to provide some guarantees in terms of determinism and robustness in order to ensure the proper operation.The adoption of mmWave frequencies for these applications may be a key enabling technology that provides a sufficient performance to cope with the most demanding scenarios, but, despite its potential, this possibility has not been completely analyzed yet.
This work aims to provide an insight of the propagation of mmWave signals in 60 GHz Industrial, Scientific and Medical (ISM) band in a typical factory workshop.To achieve this, a measurement campaign has been carried out first in a workshop using the double directional channel sounder presented in [11].Three different areas are characterized in total using 30 Line-of-Sight (LoS)/Non-Line-of-Sight (NLoS) positions, and the complete azimuth response of each position has been obtained measuring the channel impulse responses (CIRs) in 3 • steps with a bandwidth of 2.16 GHz.A stochastic channel model that represents the observed propagation has been elaborated later, where aspects as the received power, excess delay or departure/arrival directions have been modeled.The model has later been compared with two existing channel models in order to observe the differences between industrial and more ordinary scenarios.
The main contributions of this work are the following: 1) A measurement campaign of mmWave signal propagation in a workshop is carried out.Double directional Channel Impulse Responses (CIRs) from 30 positions that are divided in three zones have been acquired, and all the results have been submitted to the NIST's NextG channel model alliance public repository [12] (Measurements in a University Workshop with Industrial Equipment dataset).2) A stochastic channel model that represents the observed propagation is elaborated and validated, where the received power, delay and orientation of the multipath components are estimated based on the scenario and distance.This is available in the NIST's NextG channel model alliance public repository as well.3) An initial assessment of the performance that the 802.11ad standard can provide in industrial settings is carried out.To achieve this, the performance that the standard achieves over the proposed model is compared to the performance over other two models that represent generic indoor scenarios.The rest of the article is organized as follows.Section II presents a literature review of mmWave measurement campaigns in factory environments.All the information related to the measurement campaign carried out in this work, including details of the equipment and scenarios, is described in Section III.The postprocessing of the measured data is detailed in Section IV.Section V describes the elaborated stochastic channel model, where all the propagation parameters are detailed.The validation of the stochastic model is carried out in Section VI.In Section VII, a preliminary channel performance test can be found, where the elaborated industrial model is compared to models that describe more generic scenarios.Finally, Section VIII concludes this article.

II. EXISTING EHF BAND MEASUREMENT CAMPAIGNS IN INDUSTRIAL ENVIRONMENTS
This section discusses the existing mmWave measurement campaigns in industrial environments that can be found in the literature.Table 1 shows a summary of the main attributes of each work, and details regarding each individual work are provided next.
Several industrial scenarios as a mechanical room, experimental hall, particle accelerator ring and two maintenance tunnels are tested in [13].Here, the Talon AD7200 router [26] is employed as a 60 GHz signal source and the channel information is obtained using two criteria: by capturing IEEE 802.11ad packets and by down-converting and storing the signal.This information is used to estimate the throughput the 802.11ad standard can offer in each scenario, and a fitting to a multipath propagation model is provided as well.This model groups the multipath components in clusters, and the computed metrics are the number of clusters, number of components per cluster, delay spread (DS) and expected power of each cluster and component.
Propagation inside a high-precision machining workshop is characterized in [14] for the 3.7 and 28 GHz frequency bands.A sliding correlator double band channel sounder that provides a bandwidth of 2 GHz per frequency band is used for the measurements.A total of 76 positions are measured in total (38 LoS and 37 NLoS), and propagation metrics as path loss (PL) (fitted to the folating intercept (FI) and alpha beta gamma (ABG) models), Delay Spread (DS) and Transmitter (Tx)/Receiver (Rx) aximuth spread (AS) are extracted from them as well.
A sake factory and a subway station are measured in [15].The factory environment is composed of a 38.5 × 31.2 m 2 floor that contains a corridor and two rooms, and four Tx and five Rx positions are measured here in total.The station presents a total length of 152 m, from where the Tx and Rx are placed at three different locations each.A frequency modulated continuous wave sounder that operates at the 28 GHz spectrum is used to obtain the data from both environments.Some postprocessing is done to the data as well, from where metrics as the Path Loss (PL) (fitted to the Floating Intercept (FI) model), DS and receiver Azimuth Spread (AS) are extracted.
The PL from three different scenarios is characterized in [16] for the 24 GHz spectrum.These environments are a university workshop, a production hall and a storage room.The channel sounder employed during the measurements is composed of a signal generator, frequency band converters, omnidirectional antennas, and a spectrum analyser.A total of 26 tones that are spaced 10 MHz apart are generated for the sounding, adding up to a total bandwidth of 250 MHz, and several distances and antenna altitudes are tested in all the sites with this, considering both LoS and NLoS cases.The obtained data are used to perform a PL estimation, which is fitted to the close-in (CI) and FI models, and these parameters are also compared to the results of other works.
Measurements from a circle-shaped machine hall are presented in [17].A sliding correlator channel sounder with a bandwidth of 2 at 28 GHz is employed to obtain the CIRs.Four configurations are measured in total: two different heights and swapping the Tx/Rx ends, as one of them is static during all the process.A total of 40 positions are measured per configuration, which means that 160 measurements are provided in total.These data are processed to obtain metrics as PL (fitted to the Close-In (CI) and FI models), signal blockage, DS and AS, both in departure and arrival.
Signal propagation inside a smart warehouse filled with automatic guided vehicles is characterized in [18] for the 28 GHz spectrum.A frequency domain channel sounder is employed to carry out these measurements, which captures 1001 points in a span of 1 GHz.Three conditions are measured in total: one LoS and two NLoS, the latter being provoked by placing one and two containers.The results are processed to obtain power delay profile (PDP) and PL-related metrics.
Measurements are carried out in [19] inside a machine room with a size of 35 × 16 m 2 .A frequency domain sounder that operates in 28 GHz spectrum is employed for the measurements.This sounder uses monopole antennas in both ends, and takes 801 points spanned across 4 GHz in each acquisition.Three conditions are measured in total: none, partial, and total obstruction of the LoS component.The data are processed as well, from where metrics as PL (fitted to the FI model), DS and AS (estimated with the Space-Alternating Generalized Expectation-Maximization (SAGE) [27] algorithm) are obtained.
Propagation through an optic fiber laboratory is characterized in [20] for 28 GHz frequency band.This measurement campaign uses a pseudo-random noise (PN) channel sounder that provides a total bandwidth of 2.4 GHz, and provides two datasets.The first one measures the PL at various distances (57 positions) and the second tests several antenna orientations (12 positions).The obtained results are employed to estimate the PL (fitted to the CI and FI models), DS, and receiver AS.
In [21], sub-terahertz measurements are carried out in a corridor that separates various machining laboratories and conference rooms.CIRs are acquired via a sliding correlator channel sounder that offers a total bandwidth of 500 MHz in baseband and is centered at 142 GHz, and horn antennas are used as well to obtain double-directional responses.Thirteen measurements are carried out in total, from where 11 are LoS and 2 NLoS.These are processed later to obtain PL (fitted to the CI model), DS, and departure/arrival AS, considering the directional and omnidirectional CIRs separately as well as copolarization and cross-polarization of antennas.
Multiband measurements can be found in [22] for a production hall.This room has a total size of 171 x 74 m 2 , and four sections can be divided: three production lines and a storage area.Measurements are performed over these sections with a 15-b Pseudo-Random Noise (PN) sounder that has a total bandwidth of 6.75 GHz and is capable of acquiring the directional CIR of the 6.75, 30, and 60 GHz frequency bands simultaneously.The data are processed to estimate the PL (fitted to the Alpha Beta Gamma (ABG) model), DS and departure AS/elevation spread (ES) inside the environment.
An additional multiband measurement campaign can be found in [23].Here, propagation through a 120 × 40 m 2 textile factory is characterized using a frequency domain channel sounder.The channel sounder is configured to acquire 501 frequency points per measurement, which cover a span of 1 GHz and is centered at 3.5, 38.5, and 39.5 GHz.A total amount of 15 positions are characterized, distributed along two corridors surrounded by looms.The obtained data are used to perform PL (fitted to the ABG model), DS, and Ricean-K value estimations.
A digital fabrication facility is characterized in [24] for the 28 GHz spectrum.Here, a 13 × 13 m 2 manufacturing zone is characterized by measuring 13 positions that combine the none, partial, and total obstruction of the LoS component scenarios.The results are processed to obtain some propagation metrics as the PL (fitted to the CI model), DS, and AS of arrival.The effects of co/cross polarization of antennas are also observed for these metrics.
At last, a ray-tracing based propagation characterization is carried out in [25] for simulated light/heavy industrial sites.This work uses a simulation tool known as Wireless InSite [28] to model the propagation of signals that operate in 28 and 60 GHz mmWave spectrum using the shooting bouncing ray (SBR) technique.The outputted results include the PL, LoS probability, Rician-K factor, DS, and AS/Elevation Spread (ES), both in departure and arrival.
As it can be seen, most of the presented works are centered at 28 GHz frequency band and there are very few contributions focused on 60 GHz one.The latter is of special interest for communication systems, since the protocols 802.11ad/ay operate over these frequencies and future standards will likely use it as well.Since industrial settings contain fundamental differences compared to other generic scenarios as laboratory/office (considerably bigger rooms, dense scenarios, lots of reflectors/scatterers, etc.), existing generic measurements and models are not suitable to represent them.In this regard, several channel models as the 3GPP, NYUSIM, METIS or MiWEBA already exist in the literature, but most of them do not represent industrial scenarios.There are only two models that have recently included an Indoor Factory (InF) scenario: 3GPP and NYUSIM.However, these present some limitations, since none of them has explored 60 GHz frequency band (3GPP used measurements up to 28 GHz and NYUSIM based its metrics in a 142 GHz measurement campaign introduced in [21]), and a bandwidth as large as 2.16 GHz cannot be properly represented with them.The lack of a specific measurement based industrial channel model in 60 GHz spectrum, together with the scarcity of empirical data related to this scenario are the main motivations behind this article.

III. MEASUREMENT CAMPAIGN
This section describes the channel sounding equipment and configuration as well as the selected sites to perform the measurements.

A. CHANNEL SOUNDER
A general picture of the channel sounder is provided in Fig. 1 and all the details of its elaboration and operation can be found in [11].
This device is capable of sounding any frequency in 57-65 GHz range, with a carrier frequency step of 540 MHz.Regarding the direction, both Tx/Rx ends are equipped with  highly directional lens horn antennas that have a high frontto-back ratio, whose radiation pattern can be seen in Fig. 2.These can be mechanically steered in any orientation in the azimuth in 0.1 • steps, allowing to acquire highly directional CIRs.As for the operation, this channel sounder generates a complex baseband chirp signal that covers a total bandwidth of 1.08 GHz and is simultaneously sent to two destinations: a front-end that up-converts to the mmWave spectrum and transmits through the air interface and an oscilloscope that is used for the acquisition.The receiver front-end performs the down-conversion and forwards the signal to the same oscilloscope, and the CIR is obtained by de-convolving both Tx/Rx sequences.
In this measurement campaign, the channel sounder was configured to perform complete azimuth scans in 3 • steps for two contiguous frequency bands in the mmWave spectrum: 59.94 GHz (59.4-60.48GHz) and 61.02 GHz (60.48-61.56GHz).A total amount of 28 800 CIRs were measured per position, covering a bandwidth of 2.16 GHz.The mmWave carrier frequencies were selected to cover the channel 2 of the IEEE 802.11ad/ay standards.In order to avoid excessive acquisition times and memory requirements, each capture was limited to 10 K samples in time, which proved to provide a dynamic range of around 40 dB.Each position took approximately 12 h to sound, and the generated data occupies 1.8 GB after the CIR extraction.A summary of the most relevant parameters of the sounder for this campaign can be found in Table 2.

B. SELECTED SCENARIOS
The measurement campaign was carried out in a workshop that belongs to the university of Mondragon.This scenario is very similar to the industrial workshops presented in [14], [16], [17], where narrow corridors separate diverse sets of machines in dense environments.Three different zones were selected to perform the measurements: one that contains three vertical machining centers, another one that has three milling stations and a corridor-like environment that is between two hydraulic presses.A total of 30 different positions were acquired in total (10 per zone: 19 LoS and 11 NLoS) during Jan-Feb 2022, and the recordings were performed at night, where the machines were turned OFFand there was no human blockage.All the measurements are available in the NIST's NextG Channel Model Alliances public repository in the Measurements in a University Workshop with Industrial Equipment dataset [12].Details regarding the individual sites are provided next.

1) VERTICAL MACHINING CENTER
The vertical machining center (VMC) can be seen in Fig. 3.The environment is composed of three stations along some tool holders and metallic pillars.The wall next to the Tx node contains some metallic lockers as well.The stations have the shape of big empty boxes that are mainly composed of metallic surfaces.These elements contain doors or movable panels to access the interior, which are typically made of transparent materials as plastic.During the measurements, all these stations were empty and opened, and the tool holders remained in the same position.As for the selected positions, both LoS and NLoS scenarios were considered for the 1-6 and 7-10 Rx positions, respectively.
2) MILLING STATION Fig. 4 shows the milling station zone (Mill).This scenario is composed of three milling machines that are surrounded by have some drilling stations, two vertical machining centers and some lockers.The largest milling machine (which is below positions 17, 18, and 19) also contains a large metallic plate on its back to collect chips.The milling stations are almost entirely made of metallic elements and present an irregular shape.Their height is also similar to the one of the sounder, meaning that a partial obstruction of the LoS component is present in all the Rx positions.

3) HYDRAULIC PRESS
The hydraulic press area (HPress) can be found in Fig. 5.This scenario is composed of big and dense metallic machines that are expected to completely reflect the signals that incise on them.These machines are two hydraulic presses, an hydraulic press shear and an hydraulic press brake.Narrow corridors separate these machines and the building walls, which are located below the Tx in Fig. 5, are mostly covered in metallic  components as tools, plates, and similar.As for the positions, 3 LoS and 7 NLoS Rx locations were selected in total.

IV. DATA POSTPROCESSING
The process that extracts the relevant information from the measurements is detailed in this section.A summary of all the steps can be found in Fig. 6, and an explanation of each is provided next.
The first step of the postprocessing was to join the CIRs of the contiguous frequency bands to achieve a bandwidth of 2.16 GHz.This was possible because the channel remained invariant during the complete measurement step.
Each position was processed separately after obtaining all the corresponding CIRs.The maximum detected power was extracted from each CIR first, which contained the information of a specific antenna orientation combination.This information was used to locate the different multipath components, and a power mask was passed to the data to carry out such task.A threshold of −110 dB of relative power between the Rx/Tx signals was defined in this step, since it is the lowest value that can be applied without the interference of noise floor peaks.
As the spatial resolution of the measurements provides several samples within the half power band width (HPBW) region of the antennas, the same multipath component can be detected in contiguous orientations.These directions do not contain more than one multipath component, but rather artifacts that the antennas create when the signal travels through a direction that is not the boresight.In order to obtain the information of the real multipath component, the directions that contained energy were grouped in the next step using their proximity in delay and orientation as criteria.
An example of the processing carried out so far can be seen in Fig. 7 for the position 1 of the Vertical Machining Center (VMC) scenario.Here, a heatmap of the channel gain is displayed first using the maximum power of all the (360/3) 2 = 14400 CIRs, and the groups of the directions that are above the power threshold are shown next.The colors identify the different groups and the numbers are used to order them from the shortest to the largest delay.
The next step was to check the individual CIRs that were above the power threshold, since each one of them could contain more than one multipath component.However, each CIR showed just one multipath component in the time domain, which was attributed to the highly directional responses that are obtained with the employed antennas.Once the check for multiple multipath elements was completed, the relevant information was extracted from the previously introduced groups.To do so, the real path was assumed to be only in the direction that yields the highest power, and, consequently, the rest of the elements within the group were discarded.As for the strongest element, its power, delay, angle of departure (AoD) and angle of arrival (AoA) were stored for the stochastic model elaboration.

V. STOCHASTIC MODEL ELABORATION
The main structure of the elaborated stochastic channel model is introduced first in this section.Then, the parameterization of the model is presented, where formulas to implement aspects as the PL or delay are presented.At last, an existing implementation of the model is shown, where a table that contains all the required parameters to implement the model in any development tool is provided as well.

A. MODEL DESIGN
The starting point for the channel model design was the tapped delay line (TDL), since this model represents the propagation as a summation of multipath components of different gain, phase shift, and delay, which is an approach that is well aligned with the available data.The structure of the Tapped Delay Line (TDL) can be seen in ( 1), where β n represents the amplitude of the tap n, θ n is its corresponding phase shift and τ n is the delay β n e jθ n δ (t − τ n ) . (1) The measurements contain the directional response in the azimuth as well, which is why an extension of this model is proposed in this article.The extended TDL model can be seen in ( 2), ( where the parameters φ AoD(n) and φ AoA(n) have been included to represent the angles of departure and arrival in the azimuth plane.
In order to implement the model proposed in (2), the amount of multipath components, as well as their power, delay, and orientation must be previously known.These metrics are statistically inferred from the available measurements where formulas to represent these aspects in a stochastic or deterministic manner are elaborated.
The PL of the multipath components was observed first in individual positions, from where a great variation of the received power was detected in comparison to the free space path loss (FSPL), which is computed with the Friis transmission (3).Here, P T , P R and G T , G R represent the power and maximum antenna gains in Tx/Rx in the linear domain, respectively, λ is the wavelength of the signal in meters, and d is the physical separation between both ends in meters.
(3) This variation could be detected because the channel sounder has a dynamic range of 40 dB, and the main reason for such a wide power spread was the composition of the different surfaces that form the measured environments.The dimensions, shape, material, and orientation of the different objects that surrounded the antennas during the measurements were very diverse, where some of them reflected almost all the power to the receiver side and others acted more as scatterers, spreading the impinged signal across all the directions.
Typical channel models only consider two cases: LoS and NLoS.However, as these measurements have multipath components that can vary almost 40 dB, it was concluded that modeling all the NLoS components together would introduce a high error, and splitting this data into smaller regions was proposed as an alternative.Here, the energy arriving from all the multipath components (E s ) was compared to a Free Space Path Loss (FSPL) reference that was generated for their corresponding delay, and one of the following four categories was assigned to them depending on the result: (FSPL − 30) dB As it can be seen, four regions were defined in total.The first two are very similar, since a surface that reflects most of the impinged energy generates a gain that can be observed in the LoS case.In order to differentiate these, the antenna orientations were considered, where the directions that are close to the case of both antennas being aligned are considered as LoS and the rest are very strong reflectors.Two lower power categories were defined for the multipath components that arrived with less energy: strong and weak reflectors.Each of the presented categories has an exact range of 10 dB, which was chosen as an arbitrary delimiter that fits well with the data.With this, three power ranges that cover all the values above the −110 dB threshold that was previously introduce to detect the multipath components were defined in total.
The model design concluded with the categorization of multipath in NLoS.The next step was to elaborate the formulas that compute the amount of multipath components and their corresponding power, delay and orientation, which is discussed in the next section.

B. MODEL PARAMETERIZATION
This section describes the process carried out to generate the formulas that compute the parameters that are introduced to the extended TDL model.All the estimations were carried out for each scenario independently, where the measured LoS and NLoS positions were split as well to generate independent models for these cases.A set of parameters was estimated for the LoS path (when it was present) and the three NLoS categories: very strong, strong and weak reflectors.The process to make these estimations is detailed next.

1) NUMBER OF MULTIPATH COMPONENTS
The amount of multipath components (referred as NReflectors in this work) showed a correlation with the relative distance between stations, as the closest positions were the ones that showed the highest amount of multipath.For this reason, a linear regression was proposed to model this behavior, where the Tx/Rx separation is the independent variable and the amount of multipath components depends on it.The formula that applies the coefficients from the regression is shown in (4), where α N is the intercept point and β N represents the distance (d in meters) dependant slope.Additionally, the parameter X σ is also included in the formula to represent the regression error as a random variation of the channel.This parameter uses the mean squared error (MSE) of the regression in a random Gaussian variable that has the N (0, MSE) distribution (4)

2) PATH LOSS
Regarding the PL (which is denoted GReflectors in the model), there are two widely accepted single frequency band models in the literature that can be implemented: CI and FI.Since these models have proved to be suitable candidates to represent the propagation in countless papers, this work has analyzed their structures and selected the one that adapts best to the available data.

VOLUME 4, 2023
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.The CI model, which is also known as the log-distance model, is shown in (5).Here, FSPL( f , d 0 ) is the free space PL given by the Friis propagation (3) for any frequency f in Hz and a reference distance d 0 (typically 1 m), d represents the distance in meters, n G is the PL exponent (its value is 2 in vacuum) and X CI σ is a random Gaussian variable that represents the large-scale PL fluctuation The FI model can be found in (6), where α G is the floating intercept point in dB, β G represents the slope and, as with the CI model, X FI σ indicates the large-scale PL fluctuation with a random Gaussian variable In order to check for the suitability of the models, both the root mean square error (RMSE) and R-Squared metrics were computed for all the multipath component types in each scenario.The Root Mean Square Error (RMSE) provides the error between estimated and observed data, and it is computed as described in (7), where x i and xi are the observed and estimated samples at index i and N represents the total amount of samples.On the other hand, the R-Squared metric, which is also known as the coefficient of determination, represents how accurate are the predictions of the regression.This metric computes the proportion of the variation in the dependent variable that the regression can predict, and the result is provided as a percentage in the [0..1] range.The formula to compute this can be found in (8), where x i and xi correspond to the observed and estimated samples at index i and x is the mean.

RMSE
The performance metrics for all the scenarios can be found in Tables 3 and 4. As it can be seen, both models show similar results in the LoS path.However, this is not the case for the reflections, where the FI clearly outperforms the CI model.This is because the CI model uses the FSPL as its intercept point, which is at least 10 dB away from the weakest multipath components because this metric was used to classify them.Since the FI model is well suited to represent these cases and there is no a clear difference in the LoS case, the FI model was selected to represent the PL.

3) EXCESS DELAY ESTIMATION
As for the arrival times of the multipath components, these were calculated as random offsets of the delay that corresponds to the LoS path, and the model addresses this metric as TReflectors.The excess delays of all the detected multipath components were computed first, and the Kolmogorov-Smirnov test [30] was used in order to determine the probabilistic distribution that represents them best.Here, 24 different distributions were tested in total, from where a different best approach was proposed depending on the scenario.Even if the best approach was changing, the generalized extreme value (GEV) distribution showed a significant performance in all the tests, where its p-value was less than 0.2 below the best distribution in the worst cases.As using a different distribution for each scenario would introduce an unnecessary complexity to the model, it was decided to model all the multipath component delays using the Generalized Extreme Value (GEV) probability distribution.
The probability density function (PDF) of this distribution is shown in (9), where the three parameters that are required to define it can be differentiated: the location (μ); scale (σ ); and shape (k).An individual set of values was computed per reflector type in each scenario, as the excess delay was significantly lower in the weakest reflectors and NLoS scenarios.This distribution can generate negative values, which is not possible when generating an excess delay.As the value must be positive, the model was configured to generate numbers until a value that is greater than zero is obtained.

4) ORIENTATION IN AZIMUTH
The modeling process of the azimuth orientation, which is referred as AReflectors in the model, began by estimating the paths that each multipath component travelled through during the measurements.The layouts of the sites, as well as the Angle of Departure (AoD), Angle of Arrival (AoA), and delay of each component were used to perform these estimations, from where most components showed to be first-order reflections (i.e., signals that reflect just once before reaching the destination).Additionally, these multipath components were in angle combinations that showed a clear pattern if they were drawn in the heatmap that was presented in Fig. 7.This pattern can be seen in Fig. 8 as a blue dashed line [generated with the formula provided in (10)] in the frame that is employed to display the heatmaps ) As it can be seen, three types of regions can be differentiated in the frame: LoS, first order and second order.
The LoS region is reserved for the LoS component, and it is always centered at the (0,0) direction, which represents that both antenna boresights are aligned.
The majority of multipath components arrived in the firstorder region, which contains the directions where the signal can arrive to the destination with just one reflection.The data extracted from the measurements showed that approximately 8 out of 10 multipath components fall within this region and that they have the shortest delays, since they correspond to the shortest paths.Considering these points, the first 80% of the multipath components that the model generates are assigned to random locations within this region, and the previously introduced blue dashed line is employed as a reference in this process.The measured orientations were assumed to belong to this line, and the distance between them and the closest point of the line was used to estimate the Mean Squared Error (MSE) of the approximation.The model also generates random directions that are uniformly distributed along the line, and adds some variation to them with a random number that follows the N (0, MSE) distribution.
At last, the remaining 20% of the multipath components are assigned to the second-order region.These elements correspond to paths that contain at least two reflections, and their behavior does not follow any clear pattern as in the previous case.Since all directions seem to be equally likely, this region is modeled following a uniform random distribution.

C. IMPLEMENTATION OF THE MODEL
The previously designed model has been implemented in MATLAB under the millimeter Industrial Channel Model (mmICM) name.It is available in the NIST's NextG channel model alliance public repository under the Measurements in a University Workshop with Industrial Equipment dataset [12].The complete list of required parameters to build the model are provided in Table 5, where the values related to the number of multipath reflector types (NReflectors), their corresponding gains (GReflectors), arrival times (TReflectors), and direction in the azimuth (AReflectors) are defined for all the scenarios.This information, along the previously provided formulas, is sufficient to implement the channel model in any development tool.
The model can be configured to operate under different conditions by specifying the values of 10 parameters 1) Scenario: "VMC," "Mill" or "HPress" (default: "VMC").2) Distance: Tx/Rx separation (m) (default: 5).3) LoSPresent: "True" or "False" (default: "True").4) SamplingRate: Baseband sampling rate (SPS) (default: 2.16E9).5) SpatialRes: Spatial resolution ( • ) (default: 1).6) NSamples: Number of samples in time (default: 1E3).7) RandNumGen: "default" or any integer (default: "default").8) MinDetPow: Minimum detectable power (dB) (default: -120).9) TxAntennaRP: Radiation pattern in Tx (dBi) Coefficients must follow the spatial separation defined in SpatialRes.An empty array equals to the omnidirectional antenna (default: []).10) RxAntennaRP: Radiation pattern in Rx (dBi) Coefficients must follow the spatial separation defined in SpatialRes.An empty array equals to the omnidirectional antenna (default: []).These parameters can be specified when the channel model object is created or when a realization of a channel is Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.performed.There are four functions that allow to interact with the channel model object: 1) generateChannel(params): Performs a realization of the channel with the configuration specified in "params." The realization is divided in two main steps.The information that describes the reflectors is generated first (NReflectors, GReflectors, TReflectors, AoD, and AoA).The second step performs a convolution between each reflector and the antenna radiation pattern in order to determine the reflector gain in the spatial domain.Additional information related to the channel model can be found by introducing the help mmICM command in MAT-LAB.This command provides an example model run as well, where the same channel realization is generated with both directional/omni antennas.Fig. 9 shows both the heatmap and the CIRs that are achieved with this test.

VI. MODEL VALIDATION
The validation presented in this section has the aim of demonstrating that the channel realizations generated with the model are in accordance with the observed data.In order to achieve this, a comparison between the measurements and model generated data has been carried out as described in [31], where the cumulative distribution function (CDF) of the excess delay and dispersion in azimuth is used as the reference metric.

A. EXCESS DELAY CDF
The validation of the dispersion in time has been carried out by comparing the excess delays relative to the LoS path.This comparison has been generated for each reflector type in all the scenarios, separating the LoS/NLoS cases in two different categories.An example of such comparison for the VMC LoS scenario can be found in Fig. 10.
The measured excess delay has been compared with the averaged value of 100 realizations of the model for each scenario.As it can be seen in Table 6, the mean excess delay is slightly higher for the model in all the cases.The difference does not exceed 10% of the mean value from the measurements in most cases, with the exception of VMC LoS and Hydraulic Press Area (HPress) NLoS, where this difference is of 30% for some reflectors.These particular cases have multipath components that arrive in a wider time interval, which means that the estimations have a larger standard deviation and diverge more from the original measurements.In any case, the difference of the excess delay is never larger than 30% of the measured values which, since it is a lower deviation than the one reported in [31], it is considered to be sufficiently accurate to represent the dispersion in time.

B. AZIMUTH DISPERSION CDF
As explained in the previous section, the azimuth angles of the reflectors were selected based on the regions defined in Fig. 8.The first-order region used the curves from the figure [described in (10)] for multipath component placement, whereas the second-order uses a uniformly distributed random location from all the possible angle combinations.As the first-order region was the only one estimated with the measurements, the validation is carried out to this part alone.
The curves defined in (10) were selected for the validation of the azimuth.Here, the relative distance between the real positions and the curves was used to compute the Cumulative Distribution Function (CDF) of the dispersion.An example of this comparison can be found in Fig. 11 for the VMC LoS scenario.
The mean dispersion of the measured values was computed for all the LoS/NLoS scenarios as well as the average of 100 realizations of the model for the same conditions, the obtained results can be found in Table 7.A slightly lower dispersion was observed for the model generated data in all the cases, where the difference between the means was of 8 • in the worst case scenario.As the approach of the curves defined by (10) was based solely on first order reflections and higher order reflections do exist in the first-order region, it was concluded that this difference was reasonable enough to validate the orientation in azimuth of the model.

VII. CHANNEL PERFORMANCE ASSESSMENT
The channel performance testing aims to provide a preliminary assessment of the reliability that a commercial mmWave standard can achieve in an industrial site using the model developed in this work.To do so, the IEEE 802.11ad standard has been selected, since it is an extensively tested version and simulation tools, related literature works and commercial equipment are already available.The newer 802.11ay version was considered as well, however two main drawbacks were identified in the process.On one hand, the new mechanisms were mainly oriented to improve just the throughput, which is not a particularly demanding feature for industrial communications.On the other hand, no simulation tool or equipment was available to work with this version due to its novelty.Considering these points, the 802.11ay version was discarded and its predecessor was selected for the channel performance assessment.
This test has been carried out using the Monte Carlo simulations, where an IEEE 802.11ad physical service data unit (PSDU) is modulated and transmitted using the same channel realization with different levels of additive white Gaussian white noise (AWGN).The bit error rate (BER) obtained after performing the demodulation is stored and averaged at the end of the test, from where the Bit Error Rate (BER)/signal-tonoise ratio (SNR) curves are achieved.The WLAN toolbox from MATLAB [32] was used to perform the modem operations, as it includes functions to work with the IEEE 802.11ad standard.These simulations were carried out with the mmICM proposed in this article, as well as two additional indoor channel models that represent the mmWave spectrum, which are the 802.11ay [33] and the NYUSIM [34] models.
The 802.11ay channel model [33] follows a quasideterministic (QD) approach with ray tracing.Three different scenarios can be selected for the realizations, which are an open area hotspot, street canyon and a large hotel lobby.The 3-D environment is created first with the given dimensions, and a Tx/Rx pair is placed in a specific location.Other parameters are also used to specify antenna dimensions, beamforming, and similar settings.With this, the model would be ready to process any input signal, and this process is done by generating both deterministic (depending on the layout) and stochastic (random events) rays that are later combined to elaborate a CIR.In this work, the hotel lobby scenario was selected with a physical separation of 5 and 10 m between Tx/Rx.A single patch antenna was set in each end in order to get a radiation pattern as isotropic as possible, and both ends were faced toward each other.
The NYUSIM model [34] follows a stochastic approach that groups the multipath components in time clusters that possess different spatial lobe sizes along the CI model for PL estimation.The available scenarios for testing are urban microcell, urban macrocell, rural macrocell, indoor hotspot, and InF.There are a total of 28 configuration parameters in the model, where aspects as Tx/Rx separation, LoS presence or blockage can be defined.The model estimates the channel parameters by performing N different simulations where the receiver is placed at random locations between a predefined set of ranges.Some of the generated outputs are files that contain data related to AoD/AoA power spectrum or the Power Delay Profile (PDP).The Monte Carlo simulations carried out in this comparison are done with the configuration that approaches best tot the other models, which is the InF scenario with a fixed physical separation of 10 m and a single patch antenna in both ends.Unfortunately, the minimum distance this model operates over is of 10 m, which is why it is not possible to include its results in the 5 m simulations.
Regarding the Physical Service Data Unit (PSDU) structure for the simulation, the study presented in [10] was taken as the starting point.This work summarizes the requirements for industrial communication use cases that can be found in several studies, along with its own proposal.The most demanding scenario for Class 1 communications (which contemplates closed loop regulatory control applications) has been considered for the packet generation, where the requirements to fulfill are an update period of 0.25 ms, network size of 30 nodes and a payload of 64 B per packet.A total throughput of 30.72 Mbps was calculated for this scenario, and after analyzing the 802.11 standard [4], it was concluded that the most suitable modualtion and coding scheme (MCS) was the one, which provides a throughput of 385 Mbps along being the most robust alternative the standard can offer.The Monte Carlo testing was configured to perform 1E5 iterations in total, and Signal-to-Noise Ratios (SNRs) in the (-10:2:10) dB range were tested in each.
The obtained results can be found in Fig. 12.As it can be seen, the standard achieves a similar performance over the mmICM and 802.11ay model, at least in the 5 m simulation.This is in accordance with the measurements from both models, since a predominant LoS component is always present in the CIRs and first-order reflections are at least 10 dB lower.On the other hand, the difference that both approaches present at 10 m could be explained with the techniques that the models use to generate multipath.In mmICM, the amount of multipath components decreases as distance arises, whereas the ray-tracing approach used in 802.11ay does not consider the separation between stations and has a fixed amount of multipath.The 802.11ay model could be suffering a greater multipath at 10 m, which is why its BER is slightly higher.At last, a significantly worse performance can be appreciated with the NYUSIM channel.In this case, the measured power difference between the LoS path and reflected rays in [21] is similar to the one of the other works, but the model uses another approach to estimate the channel.The NYUSIM models the arrival of a multipath component as a group of rays that have very similar amplitudes and delays.This applies to the LoS path as well, which means that the multipath suffered in this model is significantly worse than the other ones.If this is added to the fact that the modem lacks to ability to cope with heavy multipath, the observed noise floor could be explained.At last, the error rate seen with the mmICM is similar to the one of the more generic indoor scenarios, which could indicate that mmWave communications can have the potential to operate in industrial settings.

VIII. CONCLUSION
This article proposes an industrial 60 GHz mmWave stochastic channel model elaborated from a measurement campaign, which models a setup that is mostly unexplored in the literature and provides all the function to implement it in any simulation tool.The measurements were carried out in a university workshop that was divided in three different zones.These zones contained vertical machining centers, milling stations and hydraulic presses, and 30 positions (19 LoS and 11 NLoS) were measured in total.A post processing of the measured data was performed with the aim of obtaining propagation metrics as PL, number of multipath components and dispersion in time and azimuth.A channel model named mmICM was created based on these metrics, and a validation process was carried out to verify its similarity compared to the measured values, where it was concluded that the overall behavior of the model was suitable to represent the observed environments.At last, a preliminary performance assessment of the 802.11ad standard was conducted over the model proposed in this work and two generic indoor models that can be found in the literature.Here, the expected BER was computed for several SNRs, from where a similar performance was observed in all the environments.Future works will be focused on making more advanced tests with the model, as well as expanding it with additional environments and configurations.

FIGURE 1 .
FIGURE 1.General picture of sounding system.

FIGURE 3 .
FIGURE 3. Vertical machining center zone.(a) Picture of vertical machining center.(b) Layout of vertical machining center.

FIGURE 4 .
FIGURE 4. Milling station zone.(a) Picture of milling station.(b) Layout of milling station.

FIGURE 5 .
FIGURE 5. Hydraulic press zone.(a) Picture of hydraulic press.(b) Layout of hydraulic press.

FIGURE 9 .
FIGURE 9. Example CIR generated with the model.(a) Heatmap of maximum gain per direction.(b) Boresight CIR of directional/omni antennas.

2 )
cir = generateCIR(aod,aoa,plot): Generates the complex baseband CIR of a given antenna orientation specified by the AoD and AoA parameters.The parameter plot is a boolean that generates a plot of the CIR in dB when set to "True." 3) displayHeatmap(): Displays a heatmap where the maximum amplitude of each antenna orientation combination can be seen in dB.4) displayConfig(): Displays a list with the values of the configuration parameters the model is working with.

FIGURE 12 .
FIGURE 12. Performance of an 802.11ad link evaluated with the mmICM, 802.11ay and NYUSIM models.(a) Performance at 5 meters.(b) Performance at 10 meters.