End-to-End Channel Modeling and Generation Software Based on Statistical mmWave Channel Model

Millimeter wave (mmWave) multi-path channel models based on a Saleh-Valenzuela (SV) model are fundamental for simple link-level simulations in mmWave communication systems. However, parameter extractions of the model incur developer-dependent and sometimes manual procedures, which hinders from not only confirming the validity of the existing parameters but also correcting cumulative knowledge about parameters for versatile environments. This study develops an end-to-end software to automatically perform channel modeling and simulation with statistical validity just by inputting measured power profiles. The developed software involves comprehensive algorithms to reach channel parameter extraction, such as multi-path component (MPC) extraction, MPC clustering, and maximum-likelihood estimations for respective parameters; thereby serving as a tutorial to reach the objective. Based on the developed software, we highlight that channel parameters reported by the IEEE 802.15.3c task group (TG) do not necessarily fit the original power profiles, highlighting the need to examine the validity of the existing parameters by leveraging the developed software. Moreover, we propose improved parameters for residential environments, which are reported in the TG, and demonstrate the feasibility of the developed software generating more resemble channel characteristics to the measured channel response relative to the reported parameters by the IEEE 802.15.3c TG.

mmWave communication systems. Hereinafter, we call this procedure of generating artificial CIRs based on statistical channel models as "channel generation" or "CIR generation"

A. BRIEF SURVEY OF MMWAVE CHANNEL MODELS AND MOTIVATION
Before delving into the motivation of this study, this section details the basics of the mmWave CIR generation and existing approaches that have been proposed in the literature. In more detail for the CIR generation in the mmWave band, the key channel parameters are summarized as follows [5]: (1) cluster delays, (2) in-cluster delays, (3) cluster powers, (4) in-cluster subpath powers, (5) cluster azimuth AoA, and (6) in-cluster subpath azimuth AoAs. Note that as for the fifth and sixth angular parameters, an azimuth AoD and zenith AoA/AoD should be added depending on the use cases considered. Particularly, the channel models for cellular networks require the consideration of the four angular parameters. Moreover, the number of clusters and the number of subpaths within clusters can also be modeled statistically, and in this case, the two parameters are added as the seventh and eighth channel parameters. However, to focus on the simple fundamentals, we can summarize the important channel parameters as the above six channel parameters, which are necessary to generate regardless of the use cases. Hereinafter, we refer to the above six key channel parameters as the CIR component to distinguish them from the statistical parameters of the distributions that model the CIR components.
Many approaches are proposed by research institutions and standardization communities worldwide [5], among which we focus on the mmWave channel models developed by the standardization communities or the institutions capable of large-scale measurements. First, the IEEE 802.15.3c task group (TG) contributed to developing the mmWave statistical channel models [6], [7], [8], [9] by extending the Saleh-Valenzuela model [10], which allows us to generate the above six CIR components. This standardization developed specifications for wireless personal area networks, where devices communicate with each other in close proximity, and double-directional antennas may not be required. Hence, the TG only considered the azimuth AoA. Subsequently, the IEEE 802.11ad TG extended the channel models to fit wireless local area networks [11]. Therein, the channel model included the azimuth AoD and the zenith AoA/AoD. Moreover, to statistically model the distribution of the cluster delays and cluster AoAs/AoDs, the TG applied a piecewise linear model found by ray-tracing. The SV model was applied to the other CIR components. The subsequent IEEE 802.11ay TG recommended applying ray-tracing to generate cluster delays, cluster powers, and cluster AoAs/AoDs without assuming statistical distributions [12]. The generated rays are referred to as D-ray, which are added by the residual in-cluster subpath rays, referred to as R-rays, generated by an SV model.
Outside the IEEE 802 TGs, mmWave channel models are developed, particularly for cellular networks. The 3GPP developed a statistical and map-based channel model in the technical report 38.901 [13]. These 3GPP channel models did not apply the SV models while statistical distributions on the above six CIR components including AoDs are developed so that the generated CIRs retain a pre-defined delay spread and angular spread. The mmMAGIC model from the European Union [14] and the Quadriga model from Fraunhofer HHI also applied this approach [15]. New York University (NYU) has developed a unique approach, where the two types of clusters are defined in time and space domains [16], [17]. The former cluster is called the time-cluster and is used to define the cluster delays, in-cluster subpath delays, cluster powers, and in-cluster subpath powers, while the letter cluster is called the spatial robe and is used to model AoAs/AoDs. Based on this concept, the distributions of the six CIR components are provided. For a detailed survey, readers are encouraged to refer to [5].
However, regardless of the advances in state-of-the-art channel models, the validity of the reported parameters and generated CIRs remains questionable owing to the lack of a complete algorithmic framework to extract each statistical parameter to model the CIR components and generate CIRs. Indeed, several communities provide partial procedures to extract statistical parameters (e.g., [17] from NYU); however, to the best of our knowledge, a complete framework that automatically processes the measured power delay profiles to obtain statistical parameters and generate CIRs has not been provided for any channel model.
Hence, we address this issue by developing an end-toend software to automatically perform channel modeling and generation by simply inputting measured power delay profiles by revisiting the SV mmWave channel model as an example. The SV model is chosen because the IEEE 802.15.3c TG, which also applied the SV model, provided both open datasets of measured power delay profiles and statistical parameters to model the above CIR components. This allows the verification of the developed software, which contributed to drawing the insight that the channel parameters reported in the standardizations are not necessarily fit original power profiles (see the second contribution in Section I-C). Although the developed software is specific to the SV model, this work serves as a foundation to shed light on the above issue inherent in the mmWave channel models, which is the position of this work from the viewpoint of all mmWave channel models in the literature. In the next section, we detail the previous works of the SV mmWave channel model and discuss the related issues.

B. PREVIOUS WORKS OF SV MMWAVE CHANNEL MODELS AND PROBLEMS
The Saleh-Valenzuela (SV) model offers simple and lightweight channel generations [10]. In the SV model, the inter-and intra-cluster delays and angles of arrival (AoAs) of each ray that arrives under non-line-of-sight (NLoS) conditions are modeled using canonical statistical models (e.g., Poisson arrival models for delays and uniform distributions for angles), and the powers of each ray are modeled by a simple exponential decay model. Moreover, as developed by the IEEE 802.15.3c task group (TG) as a triple SV (T-SV) model [6], [7], the SV model can be simply extended to simulate the LoS properties by adding an LoS component with the powers calculated by the two-path ground reflection model. Thus, the SV model offers basic and simple channel generations and serves as the foundation for the development of mmWave communication systems.
As for the SV model tailored for mmWave communications, several organizations have reported complete parameters. For instance, the IEEE 802.15.3c TG holistically summarized the parameters necessary for channel generation in the 60 GHz band for use in residential, office, personal desktop, and kiosk environments [8], [9]. Moreover, based on the reported parameters, TG developed the first open-source 60 GHz channel model simulator for link-level simulations with MATLAB [18], based on which, the stack of artificial CIRs is made available for each respective environment as "golden set" [19] in their terminology. Similarly, the parameters of the SV model at the 60 GHz and 80 GHz bands were reported for office and conference room environments by Laboratoire d'électronique des technologies de l'information (LETI) in France [14]. Recently, LETI further reported parameters in the frequency range of 126-156 GHz in [20]. For 60 GHz band, the IEEE 802.11ad TG also reported parameters with a particular focus on intra-cluster rays [11] based on the measurement in [21], while the intercluster characteristics are generated deterministically, not statistically, by ray tracing. The works in [22], [23], [24] also reported the parameters of intra-cluster rays.
However, the validity of the reported parameters and generated artificial CIRs remains questionable for the following two reasons. First, as identified in Section I-A, the aforementioned works only report the channel parameters and do not completely detail how these parameters are created, which hinders the validity of the reported parameters. These studies mainly focus on what is extracted, rather than how these parameters are extracted from the measured profiles. Several studies detail partial procedures to extract parameters, such as how each delayed impulse was extracted from the measured power profiles [25] or how the clustering of each impulse was performed [26], [27]. Additionally, several studies note that the maximum likelihood estimation is applied to extract parameters that well represent the characteristics of measured power profiles [28]. However, the channel parameters are final deliverables made by a process chain as an input of the measured power profile; hence, not only partial procedures but also all procedures that lie between the measurement and extracted parameters should be detailed and made programmable in public. To the best of our knowledge, there are no studies detailing the entire procedure for extracting channel parameters and performing channel generation.
The second reason is that mmWave propagation characteristics are inherently sensitive to environments, where the parameter of even one type of environment (e.g., residential, office, or conference room) cannot be represented by one concrete value. For instance, owing to the quasi-optical characteristics of mmWaves, the cluster arrival rate depends on room size. Moreover, the characteristics of the reflective paths are sensitive to the smoothness of the reflective materials [29]; hence, the power decay factors of in-cluster rays are also sensitive to these materials. For this reason, the reported parameters in the aforementioned studies may be significantly scarce, considering that many types of propagation environments exist in our living environment. Hence, a more extensive investigation and cumulative knowledge of channel parameters are necessary to conduct holistic channel simulations. To this end, programmable software that easily extracts the channel parameters given the measured power profiles is required.

C. SIGNIFICANCE, CONTRIBUTIONS, AND PAPER ORGANIZATIONS
Given this background, based on the SV model, this paper develops end-to-end software to perform channel modeling and generation automatically by inputting measured power profiles (see Fig. 2 in Section III for overview). This proposal solves the first issue by allowing the investigation of the validity of the existing parameters because all procedures to reach channel generation are clarified in a publicly programmable manner. As an example evaluation, we applied the developed software to the open-source power profile dataset provided by the IEEE 802.15.3c TG and highlighted that the parameters reported by the IEEE 802.15.3c TG do not fit the measured power profile. Based on this insight, we modified the parameters with statistical verification. Moreover, as a solution to the second problem, the developed software can be implemented on programmable hardware, such as propagation measurement systems or even communication devices. This allows us to easily collect more knowledge about channel parameters for diverse environments and generate a stack of artificial CIRs representing versatile channel environments, thereby enabling holistic channel simulations.
The salient contributions in this study are summarized as follows: • We develop an end-to-end software to perform data processing, channel modeling, and channel generation from an input of measured power profiles based on SV multi-path mmWave channel models. This proposal serves as a tutorial to guide developers smoothly in channel parameter extraction and channel generation when power profiles are available. Moreover, background problem formulations solved by each parameter extraction and their closed-form solutions are summarized to obtain the optimal parameters in the developed software. To the best of our knowledge, there are no works that provide a holistic view focusing on how mmWave channel modeling and generation can be performed, not only on what is extracted by measured channel profiles given an mmWave multi-path channel model. • As an example evaluation, we highlight that the IEEE 802.15.3c TG channel parameters do not necessarily fit the original power profiles, which gives us the insight that our developed software is required to build an accurate knowledge of the parameters. Moreover, based on the developed software, we propose improved parameters for the residential environment identified in the TG with statistical validity, and demonstrate the feasibility of the developed software to generate more accurate channel responses than the parameters reported by the TG. The remainder of this paper is organized as follows. In Section II, we review the SV model extended to mmWave communications. In Section III, we describe the algorithms used in the developed software in detail. In Section IV, we provide a numerical evaluation of our second contribution. Finally, in Section V, we provide concluding remarks.

II. REVIEW OF SV MMWAVE MULTI-PATH CHANNEL MODEL
In this section, we detail the SV model for mmWave communication. In contrast to the original SV model that considers only NLoS channel responses expressed by cluster delay, in-cluster delay, cluster power, and in-cluster power, the following extensions are generally required [6]: 1) Each impulse should also be represented by a cluster and in-cluster AoA/AoD. 2) A LoS ray should be added. The first extension allows the channel simulation to be subject to directional antennas, which are generally applied in mmWave communications. The second extension is because mmWave communications are generally considered to be used in an LoS condition because of the high attenuation owing to blockage effects [6].
Considering these extensions, the channel impulse response in this model is expressed as: where I and J (i) denote the number of clusters and in-cluster rays in cluster i. The terms T, θ AoA , ϕ AoA , θ AoD , ϕ AoD denote the delay, azimuth AoA, zenith AoA, azimuth AoD, and zenith AoD, respectively, and the subscript i denotes the realization of each respective value in the ith cluster. Moreover, their subscripts (i, j) denote the realizations of each value in the jth subpath in the ith cluster. The amplitude of each ray is denoted by A with the subscripts (i, j), where the same principle applies, except that the subscript LoS means that the amplitude is of the LoS ray. These power, delay, and AoA/AoD values were modeled using canonical statistical models, as detailed below. Among the angular variables, we only considered the azimuth AoA because the same principles were applied to the other variables. Hereafter, we refer to the azimuth AoA as "AoA" and omit the subscripts AoA and AoD in the following discussion. For both cluster and in-cluster delays, the SV model employs the Poisson arrival model. Let and λ denote the arrival rates of clusters and in-cluster rays, respectively. Therefore, the delay interval between adjacent clusters and adjacent in-cluster rays is expressed by the following exponential distribution: For both the cluster power and the power of the in-cluster rays, the SV model employs exponential power decay. As shown in Fig. 1, the cluster power, that is, the power of the first arrival ray in each cluster, is modeled as: where is the power decay factor, is the power of the first arrival ray in the first cluster, and z ∼ N(0, σ 2 1 ) is the power dispersion in the decibel domain. Second, the powers of the in-cluster rays are modeled as: where γ is the power decay factor for the in-cluster ray power, k F is the in-cluster k-factor, and z ∼ N(0, σ 2 2 ) is the power dispersion in the decibel domain.
For the cluster and in-cluster AoAs, the SV model employs uniform and Laplace distributions, respectively. These are expressed by: where σ φ is the standard deviation parameter of the Laplace distribution.
Based on the observation of the power profile, we estimate the optimal parameters of the statistical models, that is, , λ, , , k F , σ 1 , σ 2 , σ φ , which fit the measured power profile in the most likely manner. Moreover, after estimating the parameters, we can generate samples of the components of a CIR, as shown in Eqs. (2)-(7) thereby generating one sample of CIRs, as shown in (1). The entire procedure can be performed using the developed software, as detailed in the next section. Fig. 2 shows a schematic overview of the end-to-end and programmable flows to yield channel generation from the measured power profiles. This flow consists of data processing, modeling, and channel generation. The procedures are detailed in this section.

III. DEVELOPED END-TO-END SOFTWARE FOR SV MMWAVE MULTI-PATH CHANNEL MODELING AND GENERATION
Note that the outcome of the generated channel response is generally used for PHY-layer simulations, as shown in the lowest part in Fig. 2. A general procedure is: (1) filtering the response with an antenna radiation pattern of the receiver antenna of interest, resulting in a channel response only in a delay domain, and (2) convolving a transmitted waveform with the resultant channel response. This procedure results in a simulated received waveform influenced by the generated channel.

A. DATA PROCESSING PRIOR TO CHANNEL PARAMETER EXTRACTIONS 1) MPC EXTRACTION AND POWER NORMALIZATION
The measured power profiles are firstly processed by "multipath component (MPC) extraction" to eliminate a continuous power spread in both delay and angular domains; thereby extracting each delayed impulse in the measured environment. We assume that the power profiles are measured under an NLoS condition because the SV model mainly characterizes the NLoS clusters and rays therein. This is without loss of generality because, even if the measurement is conducted under the LoS condition, we can run the following algorithm, where an LoS ray can be identified. Thereafter, the LoS ray can be eliminated.
In the power profiles, as shown in the input in Fig. 2, a continuous power spread in the delay time domain always exists because the measurement bandwidths are finite. Moreover, for the angular domain, a continuous power spread exists because the receiving antennas rotated during the measurement possess power gains not only in the pointing direction but also in other angles of arrival. This fact can be mathematically expressed as: where a k ∈ R + , φ k ∈ [0, 2π), τ k ∈ R + ∪ {0}, and θ k ∈ [0, 2π) denote the amplitude, phase, delay, and azimuth AoA of the kth impulse, respectively. Operator · takes the magnitude of a complex number. The terms H(τ )∈ C and A RX (θ )∈ C denote the joint filters of the transmitter, receiver, and antenna pattern of the receiver, respectively. Note that in A RX (θ ), θ = 0 is the antenna boresight without loss of generality, meaning that the A RX (0) 2 ≥ A RX (θ ) 2 ∀θ . Moreover, we also assume H(τ ) is centered at τ = 0 without loss of generality, meaning that H(0) 2 ≥ H(τ ) 2 ∀τ . In summary, the first step is to deconvolve the measured power profiles and estimate each delayed impulse in the measured environment.
It should be noted that phase φ k is not usually estimated if channel generation is targeted as the objective of the measurements. This is because, in the channel generation at the end, the phase is given by a sample of uniform distribution owing to the sensitivity of the phase subject to the path distances. This implies that this phase term is not modeled based on real values from the measurement; hence, we omit the estimation of phase φ k . Thus, the problem boils down to estimating the amplitude, delay, and AoA of each delayed impulse, and the triplets (a k , τ k , θ k ) K k=1 are usually called "parameters of MPC," or simply MPC. Hereafter, we use the terminology "delayed impulse" and "MPC" interchangeably, while we mainly use the latter to avoid confusion.
To this end, we employ successive interference cancellation, which is the most basic algorithm for extracting MPC. The algorithm employs the following two facts: First, the power spread in the proximity of (τ k , θ k ) in P(τ, θ ) basically forms a mountain shape with the peak of a 2 k at the coordinates of (τ k , θ k ). This can be easily proven by the fact that the power profile for each delay and angle P(τ, θ ) in Eq. (8) is calculated as follows: based on the triangle inequality. Recalling that A RX (0) 2 ≥ A RX (θ ) 2 ∀θ , and H(0) 2 ≥ H(τ ) 2 ∀τ , the above statement can be proven. Hence, by searching for the power peak in power profile, we can easily obtain one MPC. Secondly, the filter patterns and antenna patterns in the power domain can be known prior to measurements by referring to the specifications of the filters and antennas, or by conducting an anechoic room measurement. Hence, the continuous power spread can be eliminated, and the resultant power profile is without the MPCs previously found. This power profile can be re-used to search for other MPCs.
To summarize, the successive interference cancellation is an iterative algorithm, where the maximum power peak is searched, the power, delay, and angle of the peak are stored, and the continuous power spread is eliminated by the known joint power pattern of the filter and antenna gain, k . This algorithm is summarized in Algorithm 1. H re and PADP indicate the reference pattern of the continuous power spread and an observed power profile respectively. P th indicates the threshold of the MPC power and is set to be +20 dB above the noise floor, which is sufficient to ensure the standard recommendation of +18 dB in [30]. Note that this procedure is not limited to power delay profile, that is, if one measures the channel response expressed inside the operator · 2 in the right-hand-side in (8), we can perform the successive interference cancellation by using H(τ )A RX (θ ), which can be measured in an anechoic room.
Subsequently, the power of the extracted MPCs is normalized to eliminate the power dependency with respect to the transmitter-receiver distance. This is performed by subtracting a 2 k from the path loss calculated by the distance between the transmitter and receiver in which the power profile is Algorithm 1 MPC Extraction and Power Normalization input: H re , PADP, pathloss calculated by the transmitterreceiver distance in the measurement, loss output: measured. In this step, arbitrary standard path-loss models that fit scenarios of the measurement (e.g., the 3GPP path loss model in the 38.901 channel report [13]) can be applied.

2) MPC CLUSTERING
After extracting the MPCs, we perform clustering for the extracted MPCs with respect to the similarity of their characteristics in terms of the delay τ and AoA θ . A naive approach is to cluster them manually through visual inspection. However, one of our objectives is to provide a programmable flow throughout the procedure; hence, we apply a K-power means algorithm [26], which is an automatic clustering algorithm.
Briefly, the K-power means algorithm [26] searches for optimal "centroids" to cluster the MPCs, thereby generating a clear separation between distant MPCs in the delay and AoA space (τ, θ ). The centroid is defined as points in a delay and AoA space (τ, θ ), where the number of centroids is equal to the number of clusters. Once the centroids are initialized, they determine the initial cluster of MPCs, where an MPC is clustered in the ith cluster if and only if the MPC and the centroid are the closest among all initialized centroids. In addition, these centroids are iteratively moved to minimize the sum of the distances between each centroid and the MPCs that belong to the centroid. This procedure automatically yields clusters of MPCs, where two distant MPCs in the delay and angle spaces are categorized into separate clusters. Note that in the K-power means algorithm, the distance in the delay and angle space is defined as the Euclidean distance weighted by the power of each MPC, so that MPCs with larger powers tend to be categorized into separated clusters. Store current centroids c i ← c i ∀i # Cluster MPCs based on current centroids 5: for each MPC index k do 6: for each cluster index i do 7: Set MPC coordinate p k ← (ατ k , sin θ k , cos θ k ) Algorithm 2 summarizes the K-power means algorithm. The algorithm follows the original K-power means algorithm [26], with an exception that it supports angular circularity. This is simply achieved by defining the space as (ατ, sin θ, cosθ) instead of (τ, θ ), where α is a weight constant. In Algorithm 2, c 1 , . . . , c k ∈ R 3 and I denote the initial cluster centroids and the number of clusters, respectively.(i k ) K k=1 indicates the cluster indices of each MPC. · 2 indicates the L2 norm of the vector. In addition, ε is the threshold to decide the algorithm convergence and is set to 10 −6 multiplied by the element that possesses the maximum value among the initialized centroids c 1 , . . . , c k . Moreover, α = 6/(max τ k − min τ k ) denote the converging factor and the scaling factor, respectively. ℵ i indicates the MPC indices assigned for kth cluster. τ C , ψ sin , and ψ cos are the normalized cluster centroids of the delay and triangle functions of AoAs, respectively.
Note that prior to the estimation of the optimal centroids, the initial cluster centroids are determined by the max-min algorithm, which is also summarized in Algorithm 3 [31]. This is intended to enhance the characteristics of the converged centroids in Algorithm 2 such that the centroids Algorithm 3 Max-Min Algorithm input: MPCs, (a k , τ k , θ k ) K k=1 , # clusters, I output: initial cluster centroids c 1 , . . . , c I 1: Set MPC coordinate p k ← (ατ k , sin θ k , cos θ k ) 2: Select two MPC coordinates with largest distance as c 1 and c 2 3: for each cluster index i = 3, 4, . . . I do 4: for each cluster index i = 1, 2, . . . i − 1 do 5: Among (p k ) K k=1 , select MPC coordinate with minimum distance to c i as p (i) i as c i 9: end for 10: return c 1 , . . . , c I are dispersed in the space of (ατ, sin θ, cosθ); thereby preventing the MPCs close in the space from being categorized in the same cluster.

B. PARAMETER EXTRACTION ALGORITHM AND BACKGROUND PROBLEM FORMULATIONS 1) CLUSTER DELAY AND POWER PARAMETER EXTRACTION
Firstly, we estimate the parameters related to the cluster delay and powers, that is, the cluster arrival rate , the first-cluster power cluster power decay factor , and standard deviation of cluster power dispersion σ 1 . The entire algorithm is shown in Algorithm 4, which consists of a preparation step, cluster arrival rate estimation step, and a step to estimate other parameters, which are detailed as follows: 1. Data Preparation: First, we extract a representative MPC in each cluster; thereby extracting the cluster delay, cluster AoA, and cluster amplitude. This is done by extracting the MPC with the largest power for each cluster i = 1, 2, . . . I, which is formulated as: where (a (i) , τ (i) , θ (i) ) I i=1 is the cluster amplitude, cluster delay, and cluster AoA for each cluster, respectively, and are used for the following parameter estimation.

Cluster Arrival Rate Estimation:
The most-likely cluster arrival rate can be estimated based on the fact that the delay interval between successive clusters follows exponential distribution with the rate parameter of . This problem can be formulated as Extract representative MPC as in (10) and (11) 3: end for # 2. Cluster arrival rate estimation 4: Calculate the optimal cluster arrival rate as in (13) # 3 Estimation of other parameters 5: Calculate optimal ln as in (20) 6: Calculate optimal inv as in (21) 7: Calculate optimal σ 1,ln as in (22) 8: Calculate optimal first cluster power ← e ln 9: Calculate optimal power decay factor ← −(1/ inv ) 10: Calculate optimal std. of cluster power dispersion σ 1 ← (10/ ln 10)σ 1,ln 11: return , , , σ 1 Note: std.: standard deviation.
where the objective function is the likelihood given the rate parameter , that is, The optimal rate parameter can be determined by solving the optimization, which is proven to be the reciprocal of the sample mean: Hence, as an algorithm implementation guideline, we can simply calculate the above reciprocal of the sample mean to automatically extract the most-likely rate parameter. This function is implemented in Algorithm 4.

Estimation of Cluster Power Decay Factor, First-Cluster Power, and Standard Deviation of Cluster Power Dispersion:
Parameters related to the cluster power, that is, the firstcluster power cluster power decay factor , and the standard deviation of cluster power dispersion σ 1 are also estimated using the maximum likelihood method. To explain the background problem formulation for estimating the mostlikely parameters, we reformulate the logarithm of the cluster power model in Eq. (4) as: where ln = ln( ), Z ln = Z 10 ln 10 ∼ N 0, σ 2 1,ln Based on the reformulation, the problem narrows down to estimating the most-likely parameters ln , inv , and σ 1,ln given the measured samples of cluster powers and delays (a 2 (i) , τ (i) ) I i=1 . Based on these equations, we draw the problem formulation of the maximum-likelihood estimation, solution of the problem, and algorithm implementation guidelines. From Eqs. (14) and (17), each observed cluster power (ln a 2 (i) ) I i=1 should follow the normal distribution N( ln + inv (τ (l) − τ (1) ), σ 2 1,ln ). Hence, the likelihood of the observation (ln a 2 (i) ) I i=1 is derived as which should be maximized with respect to the parameters ln , inv , and σ 2 1,ln . The problem is formulated as follows: maximize ln , inv ,σ 2 1,ln (19).
By solving this problem, the most likely parameters can be extracted in the following closed form: Hence, as the algorithm implementation guideline, we can simply calculate ln , inv , and σ 1,ln from the above equations and calculate the optimal variables , , and σ 1 by inversely solving (15), (16), and (18), respectively. These calculations are implemented in Algorithm 4.
It should be noted that the above solutions are equivalent to the outcome of least-squared fitting [32]. Hence, one can simply replace the calculation part for obtaining ln , inv , and σ 1,ln with the least-squared fitting function defined in libraries for respective programming languages (e.g., "lsqr" function in MATLAB [33]). However, one of our objectives is to provide a tutorial on the methodology to extract each parameter while detailing the background problem formulation solved in this parameter extraction part. Hence, we provided the guidelines by starting the problem formulations, and subsequently, by implementing the closed-form solutions as in (20)- (22).

2) IN-CLUSTER DELAY, POWER, AND AOA PARAMETER EXTRACTION
Secondly, we estimate the parameters related to the delay, AoA, and power of in-cluster rays, that is, in-cluster arrival rate λ, in-cluster AoA dispersion term σ φ , in-cluster Algorithm 5 In-Cluster Delay, Power, and AoA Parameter Extraction input: MPCs, (a k , τ k , θ k ) K k=1 , Cluster indices for each MPC (i k ) K k=1 , representative MPCs in each cluster a (i) , τ (i) , θ (i) I i=1 , optimal cluster decay factor , optimal first-cluster power output: In-cluster arrival rate λ , in-cluster AoA dispersion term σ φ , in-cluster k-factor k F , in-cluster power decay factor γ , and std. of in-cluster power dispersion σ 2 # 1. Data Preparation 1: for each cluster index i = 1, 2, . . . I do 2: Filter MPCs as in (23) 3: Sort MPCs with the delay value as in (24) and (25) 4: Set in-cluster amplitude, delay, and AoA as in (26) 5: end for # 2. In-cluster arrival rate estimation 6: Calculate the optimal in-cluster arrival rate λ as in (28) # 3 In-cluster AoA dispersion term estimation 7: Calculate the optimal in-cluster AoA σ φ as in (30) # 4 Estimation on other parameters 8: Calculate optimal k F,inv as in (36) 9: Calculate optimal γ inv as in (37) 10: Calculate optimal σ 2,ln as in (38) 11: Calculate optimal in-cluster k factor k F ← −k F,inv 12: Calculate optimal in-power decay factor γ ← −(1/γ inv ) 13: Calculate optimal std. of cluster power dispersion σ 2 ← (10/ ln 10)σ 2,ln 14: return λ , k F , γ , σ 2 , σ φ Note: std.: standard deviation. k-factor k F , in-cluster power decay factor γ , and standard deviation of in-cluster power dispersion σ 2 . The entire algorithm is shown in Algorithm 5, which consists of a preparation step, cluster arrival rate estimation step, AoA dispersion estimation step, and a step to estimate other parameters related to the power of the in-cluster rays, which are detailed as follows: 1. Data Preparation: As data preparation, we prune the MPCs with either following property: 1) The MPCs with a lower delay value than the cluster delay τ (i) ; 2) The MPCs with a lower power than a pre-defined threshold P th2 .
Pruning MPCs with the first property ensures the consistency with the SV mmWave channel model. Recalling that the model considered one-sided exponential decay starting from the cluster arrival time τ (i) , the MPC with a delay value lower than τ (i) is a nuisance. Pruning the MPC with the second property allows an arbitrary graduality in the generation of in-cluster rays. A lower P th2 allows the generation of more in-cluster rays, which results in more fine-grained and accurate CIRs with the risk of a higher computational complexity in subsequent PHY layer simulations. Meanwhile, a higher P th2 allows the generation of limited in-cluster rays with a higher power, which results in a coarse-grained CIR and is beneficial in terms of computational complexity in subsequent PHY layer simulations.
For each cluster index i = 1, 2, . . . I, the MPC indices that satisfy the following conditions are searched for: Subsequently, the elements in K i are sorted in the order of the delay value as follows: subject to: where J (i) = |K i | is the number of MPCs in cluster i found in this step. Finally, the amplitude, delay and AoA of the MPCs used in the following parameter extraction for each cluster i = 1, 2, . . . I are obtained by: Note that the delay and AoA are standardized by the cluster delay τ (i) and cluster AoA θ (i) , respectively, for subsequent parameter extraction.
2. In-Cluster Arrival Rate Estimation: Similar to the estimation of cluster arrival rate, the most-likely in-cluster arrival rate can be estimated based on the fact that the delay interval between successive clusters follows exponential distribution with the rate parameter of λ. This problem can be formulated as a maximum-likelihood estimation problem as follows: (27) where the objective function is the likelihood of the observations τ (i,j+1) −τ (i,j) for i = 1 · · · I and j = 2, . . . J (i) given the rate parameter λ, that is, The optimal parameter can be obtained in the following closed form: Hence, the implementation guideline calculates the optimal in-cluster arrival rate by using the above equation. This is also implemented in Algorithm 5.
3. In-Cluster AoA Dispersion Estimation: The most-likely in-cluster AoA dispersion term can be estimated based on the fact that the in-cluster AoA follows an exponential distribution with the parameter of σ φ .
where the objective function is the likelihood of the observations θ (i,j) for i = 1 · · · I and j = 2, . . . J (i) given the parameter σ φ , that is, The optimal parameter is expressed in the following closed form: This calculation is implemented in Algorithm 5.

Estimation on Other Parameters Related to In-Cluster Ray Power:
The most likely parameters related to the in-cluster ray power are also estimated using the maximumlikelihood estimation. Similar to the estimation of the cluster powers, we reformulate the logarithm of the in-cluster ray power model in (5) to clarify the background problem as follows: where Note that for the first cluster power and cluster power decay factor in the model, we use the optimal values found by Algorithm 4, which appears as the input in Algorithm 5. From the above reformulation, the problem boils down to estimating the optimal values for γ inv , k F,inv , and σ 2,ln that maximize the likelihood function of the in-cluster subpath power observations ln a 2 (i,j) for i = 1 · · · I and j = 2, . . . J (i) . Based on this reformulation, we can calculate the likelihood of the in-cluster subpath power observations, based on which we can draw the expression for the optimal parameters and the implementation guideline. For the in-cluster subpath power observations ln a 2 (i,j) for i = 1 · · · I and j = 2, . . . J (i) , the likelihood function is expressed as: which should be minimized with respect to γ inv , k F,inv , and σ 2,ln . The problem is formulated as follows: .
By solving this problem, the most likely parameters can be extracted in the following closed forms: where (1) .
Hence, as the algorithm implementation guideline, we can simply calculate k F,inv , γ inv , and σ 2,ln from the above equations and calculate the optimal variables k F, , γ , and σ 2 by inversely solving (31), (32), and (34), respectively. This calculation is implemented in Algorithm 5. Note that the above solutions are equivalent to the outcome of leastsquares fitting [33] with (τ (i, j) ) i, j and (B (i,j) ) i, j . Hence, this calculation can be replaced with the least-squares fitting function for the respective programming language.

C. CHANNEL GENERATION PROCEDURE
In this final step, a realization of the CIR is generated based on the estimated parameters and the probabilistic models described in Section II. First, we set the optimal parameters found in the previous step into the probabilistic models for the CIR components, that is, cluster delay, in-cluster delay, cluster AoA, in-cluster AoA, cluster power, and in-cluster subpath power. Thereafter, each value of the CIR component is sampled. Finally, the LoS component is added in both the delay time and AoA of zero with the unit amplitude. Algorithm 6 summarizes the entire flow of the CIR generation.
Note that the power of the LoS component can be replaced by that calculated by the ground reflection model, which results in the T-SV model proposed in the IEEE 802.15.3c TG [6]. This modification yields more accurate CIRs, particularly when a desk usage scenario is considered. In this case, the amplitude of the LoS component is calculated as where R 0 is the reflection coefficient, λ w is the wavelength, h TX and h RX are the heights of the transmitter and receiver, respectively, and L is the distance between them. This replacement should occur when the direct wave and reflected wave from a desk cannot be resolved owing to a lack of measurement bandwidth. This point was identified by IEEE 802.15.3c because one of the use cases was a personal area network in a desktop scenario.

16:
Calculate in-cluster subpath power A 2 Note: Poi(·) indicates the Poisson distribution with the input of the average number of arrivals. Uniform(a, b) for a, b∈ R, a < b indicate the uniform distribution whose sample ranges from a to b. Laplace(0, c) for c∈ R indicate the Laplace distribution with the mean value of zero and the scale value of c.
The data is provided by National Institute of Information and Communications Technologies (NICT), and based on the data the parameters of the channel model were extracted in the IEEE 802.15.3c TG. This dataset is available in [34] and is detailed below.
The measurement campaign for the acquisition of the dataset was conducted in a "LoS residential" environment, with the specifications shown in Table 1 [34]. The measurement was conducted with a vector network analyzer (VNA) to obtain frequency-wise channel state information with a frequency step of 7.  following two types of antennas were employed: an omnidirectional antenna in the azimuth plane and a directional antenna with an HPBW of 15 • . Power profiles were obtained for each antenna type, and accordingly, the parameters of the channel model were recorded for each transmitter antenna.
Before the comparison with the IEEE 802.15.3c channel report, we checked the statistical validity of the parameters extracted by the developed software for the statistical models in eq. (2)-(7). This validation was performed with the MPC dataset extracted by applying Algorithm 1 to the aforementioned dataset with the transmitter HPBW of 15 • . Moreover, the MPCs were clustered with Algorithm 2. This validation result is detailed and visualized in Appendix.

B. EXTRACTED PARAMETER COMPARISON WITH IEEE 802.15.3C CHANNEL REPORT 1) PARAMETERS FOR CLUSTER DELAY AND CLUSTER POWER
We analyze the extracted parameters related to the modeling of the cluster delay and cluster power, that is, cluster arrival rate , cluster power decay factor , first-cluster first-ray power gain , and cluster power decay shadowing standard deviation σ 1 . Table 2 shows the comparison between the parameters reported in the IEEE 802.15.3c and the optimal parameters extracted by the developed software. As shown in Table 2, several parameters result in a mismatch with the optimal parameters extracted by the developed software. Recalling that both use the same dataset, and that the developed software calculates the optimal parameters to fit the data, the parameters reported in the IEEE 802.15.3c TG are not optimal for fitting the data.
To visualize the fact that the parameters reported by the IEEE 802.15.3c TG does not fit the original data, we show the power decay line and the standard deviation of the power dispersion from the line for the transmitter antenna with an HPBW of 15 • in Fig. 3. The power decay line based on the IEEE 802.15.3c channel report largely deviates from the data, particularly for the second and fourth clusters. Accordingly, the power decay line based on the IEEE 802.15.3c channel report overestimates the power range within the standard deviation depicted in the shaded area. This leads to an  inaccurate CIR generation, considering that in the channel generation step, the cluster power is basically generated inside the area within the standard deviation. By contrast, the power decay line extracted by the developed software fits well with the MPC data owing to the maximum-likelihood estimation summarized in Section III-B1. Hence, we can conclude that the parameters reported in the IEEE 802.15.3c are not optimal and can be modified by the developed software.
To analyze the cluster arrival rate , the values extracted by the developed software and those in the IEEE 802.15.3c channel report are close, and in both results, the arrival rate in Tx HPBW 15 • is smaller than that in Tx HPBW 360 • in Table 2. The intensity of the cluster arrival represents multipath richness; hence, these results are consistent with the propagation mechanism, where the narrower beam results in fewer clusters.

2) PARAMETERS FOR IN-CLUSTER DELAY, POWER, AND AOA
Second, we analyze the extracted parameters related to the modeling of the in-cluster delay, in-cluster ray power, and in-cluster AoA, that is, the in-cluster arrival rate λ, in-cluster power decay factor γ , in-cluster k factor k F , in-cluster power decay shadowing standard deviation σ 2 , and in-cluster AoA dispersion term σ φ . As stated in Section III-B2, these parameters are sensitive to the intensity of the pruned in-cluster MPCs during the data preparation step. Moreover, as is shown later, the IEEE 802.15.3c TG also pruned the MPCs to extract these parameters. Hence, we first analyze the dependency on the power threshold P th2 , thereby predicting the power threshold to reach similar parameters reported in the IEEE 802.15.3c. Subsequently, we compare the resultant parameters with those reported in the IEEE 802.15.3c TG.
Dependency on MPC Pruning Threshold P th2 : The threshold for MPC pruning P th2 is expected to have a dominant effect, particularly on in-cluster ray arrival rate and in-cluster k-factor. This is because a large threshold results in a coarse arrival in the in-cluster subpath, which lowers the in-cluster subpath arrival rate. Moreover, recalling that the in-cluster k factor is a bias term in the fitting line represented in (30), a large threshold results in a larger value of k F,inv and a smaller value of k F .
For this reason, we investigate the effect of P th2 on the extracted in-cluster arrival rate λ and in-cluster k-factor k F as shown in Fig. 4. As predicted above, a larger threshold, that is, more intensive MPC pruning with lower powers, results in a lower in-cluster ray arrival rate and in-cluster k factor. More importantly, by the IEEE 802.15.3c TG, a lower in-cluster arrival rate and in-cluster k-factor are reported, implying that the TG also pruned MPCs with a threshold as a hidden assumption. Hence, to ensure a fair comparison, we select the threshold yielding parameters close to those reported in IEEE 802.15.3c, that is, P th2 = +45 dB above the noise floor for the case of the 15 • HPBW transmitter antenna and P th2 = +29 dB above the noise floor for the case of the omnidirectional transmitter antenna above the noise floor.
Parameter Comparison: Table 3 shows the parameters related to the in-cluster properties reported by the IEEE 802.15.3c TG and extracted by the developed software. Note that for the developed software, the aforementioned MPC pruning threshold is used, and for this reason, the in-cluster subpath arrival rate λ and in-cluster k-factor extracted in the developed software are close to those reported by the IEEE 802.15.3c TG. Nonetheless, we can observe a mismatch in the other parameters. In particular, the IEEE 802.15.3c TG overestimates the in-cluster AoA dispersion term σ φ relative to the developed software, and this validity should be verified in comparison with the measured in-cluster AoA detailed next.
As an example, Fig. 5 shows the cumulative distribution function of the in-cluster AoAs in the measurement and those of the Laplace distribution with parameters in the IEEE 802.15.3c TG and in the developed software for the 15 • HPBW of the transmitter antenna. The parameters extracted by the developed software fitted well with the distribution of the measured in-cluster AoAs. Hence, we can also conclude that the parameters reported in the IEEE 802.15.3c are not optimal, and that the parameters can be modified by the developed software.

C. COMPARISON OF THE GENERATED CHANNEL RESPONSE
Finally, we evaluate the channels generated by the developed software, thereby showing the impact of the optimality of the developed software shown in Figs. 3 and 5 on the generated channels. Fig. 6 shows the measured MPC and generated MPCs in delay and power ((a)-(c)), and in delay and AoA ((d)-(f)). We also generated MPCs based on the IEEE 802.15.3c channel report by running Algorithm 6 after setting the reported channel parameters. However, exceptionally, we set the cluster delay and AoA to be the same variables as those in the measurement for both the developed software and the IEEE 802.15.3c channel report to focus on the in-cluster characteristics.
Impact on delay vs. power: From Figs. 6 (a)-(c), we can see that the MPCs generated based on the IEEE 802.15.3c channel report resulted in a larger power in the third cluster than the measured MPCs. This is attributed to the characteristics shown in Fig. 6, where the IEEE 802.15.3c channel report overestimated the power range between −σ 1 and +σ 1 relative to the fitted curve. By contrast, the MPCs generated by the developed software do not suffer from this overestimation, resulting in a more realistic channel generation in the sense that the generated MPCs are close to the measured MPCs.
Impact on delay vs. AoA: From Figs. 6 (d)-(f), we can see that the MPCs generated by the IEEE 802.15.3c channel report resulted in a much larger in-cluster AoA dispersion than the measured MPCs. This is attributed to the aforementioned property of the in-cluster AoA dispersion term σ φ shown in Fig. 6, where the distribution of the in-cluster AoA possesses a larger tail for both sides. This characteristic is far from that of the measured MPCs, as shown in Figs. 6(d) and (e). Meanwhile, the MPCs generated by the developed software did not yield such a large dispersion, leading to the similar characteristics to the measured MPCs in the in-cluster AoA.

V. CONCLUSION
We presented an end-to-end channel modeling and generation software for mmWave communication systems. The software showed programmable flows for data processing, parameter extraction, and channel generation. The software was validated using open-source mmWave measurement datasets. Consequently, the multi-path propagation characteristics can be described using a set of parameters and reproduced with the extracted parameters. Owing to its transparency and simple implementation, the developed software is expected to be applied in the system development and standardization of near-future mmWave communication systems.

APPENDIX
As explained in Section IV-A, we checked the validity of the statistical distributions with the extracted parameters for each variable in eq. (2)- (7). For this validation, we used the MPC data that was extracted from the power delay profile dataset with the transmitter HPBW of 15 • [34] by applying Algorithm 1. Moreover, the MPCs were clustered with Algorithm 2.   Fig. 7 shows the statistical distributions for each variable in eq. (2)-(7) with the empirical distributions of corresponding variables found from the MPC data. Note that the parameters are found by Algorithms 3 and 4. To apply Algorithm 4, we adopted the MPC pruning threshold P th2 = 45 dB above the noise floor. The models with the VOLUME 4, 2023 837 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. extracted parameters do not contradict the empirical distributions of the MPC data due to the maximum likelihood distributions formulated in Section III.

ACKNOWLEDGMENT
Moreover, the authors thank Dr. Hirokazu Sawada in NICT for a valuable advice for the channel modeling.