Channel Estimation for Diffusive MIMO Molecular Communications

Channel Estimation for Dif fusi v e MIMO Molecular Communications S. Mohammadreza Rouzegar , Umberto Spagnolini Dipartimento di Elettronica, Informazione e Bioingegneria, Politecnico di Milano, Milan, Italy Email: seyedmohammadreza.rouze gar@mail.polimi.it, umberto.spagnolini@polimi.it Abstract —In diffusion-based communication, as f or molecular systems, the achievable data rate is v ery low due to the slow nature of diffusion and the existence of sever e inter -symbol- interference (ISI). Multiple-input multiple-output (MIMO) tech- nique can be used to improv e the data rate. Knowledge of channel impulse response (CIR) is essential for equalization and detection in MIMO systems. This paper presents a training-based CIR estimation for diffusive MIMO (D-MIMO) channels. Maximum likelihood and least-squares estimators are derived, and the training sequences are designed to minimize the corresponding Cramér -Rao bound. Sub-optimal estimators are compared to Cramér -Rao bound to validate their perf ormance. Index T erms —Molecular communication, Diffusi ve Multiple- input Multiple-output, Channel impulse response, Cramér -Rao Bound, T raining Sequence Design I . I N T R O D U C T I O N Molecular communication (MC) is a bio-inspired solution of communication at nano-scale [1], [2]. Con ventional com- munication systems transfer information using electromagnetic wa ves. At nano-scale, antennas suffer the constraint of being at comparable scale of electromagnetic wav elength. Addi- tionally , using electromagnetic wave for nanomachines can be detrimental in some en vironments, such as inside a body where electromagnetic radiation can be harmful for health. Hence, MC can be a preferred solution for communication among nanomachines to b uild a nanonetwork, so they perform complex tasks which could not be possible individually [3], [4]. In MC, bio-nanomachines communicate through exchang- ing molecules. In fact, the simplest system needs a trans- mitter to send the information molecules, and a recei ver to collect them. Information molecules dif fuse tow ard the receiv er using Brownian motion resulting from their collision with the molecules in the fluid [5], [6]. Information can be encoded to the dif ferent properties of molecules, such as their concentration [7], number [8], type [9], and time of release [10]. One of the main challenges of MC is to deal with the long tail of diffusi ve propagation that causes sev ere inter-symbol- interference (ISI). One can increase the symbol interv al time to eliminate the ISI, but the need for higher data rate justifies the optimization of the symbol interv al time to hav e few channel taps due to the ISI [11]. Even if one optimizes the symbol interval time, the slow nature of diffusion makes the data rate still low . Using multiple-input multiple-output (MIMO) technique is a widely in vestigated solution to address this problem [12], [13] and can be adopted for MC. In this paper , we assume that information is encoded in the number of molecules observed at the receiv er . Therefore, we define the channel impulse response (CIR) as the expected number of molecules observed at the receiver domain at time t denoted as ¯ c ( t ) after instantaneous release of molecules at t = 0 [14]. W e assume the number of molecules at the receiv er follows the Poisson distribution as introduced in [14], [15] and the CIR is their mean number of molecules. The goal of this paper is to gain insight to the estimation of the CIR of D-MIMO communication channel. In [13], authors in vestigated various div ersity technique in D-MIMO communication assuming full knowledge of CIR. The authors of [12], modeled the 2 × 2 D-MIMO channel by fitting a curve to the simulated data. Goal of this paper is to define a training- based channel estimation for D-MIMO molecular communi- cation. In [14], authors introduced the channel estimation for the Poisson MC channel and proposed the estimators of CIR in single transmitter and single receiver system. Nov elty of this work can be considered the introduction of a system model and notation (Section 2) that allo w us to estimate the CIR of M × M D-MIMO Poisson channel. W e extended the steps of R. Schober et. al. [14] to D-MIMO channel estimation (Section 3) by accounting for the inter-link diffusi ve interference. Furthermore, we propose in Section 4, a D- MIMO specific method for designing the training sequence that minimize the Cramér-Rao bound (CRB) at all receivers simultaneously . I I . S Y S T E M M O D E L W e consider a M × M D-MIMO system for MC as sho wn in Fig. 1. The system consists of M pair of transmitters denoted as T x i , and receiv ers Rx j , where i, j ∈ { 1 , 2 , 3 ...M } . W e assume that all transmitters emit the same type of molecules. Each transmitter emits a known number of molecules N at the beginning of each symbol intervals. The molecules diffuse in the environment and some of them reach the M receiv ers. Transmitters and receiv ers are not fixed in their position and could slightly mov e on fluid where molecules diffuse, so the CIR changes o ver time. W e assume a block-type communication and we estimate the CIR at the beginning of ev ery block by sending a properly designed training sequence, and we assume that the D-MIMO channel does not v ary during Fig. 1. T opological model for M × M D-MIMO system. The M transmitters ( T x 1 , ..., T x M ) release the same type of molecules (pentagon shape), and molecules here have different colors according to the corresponding transmitter . each block. Then, the estimated CIR is used for equalization at the receiver ov er the rest of the block. T ransmitters modulate the molecules density using concen- tration shift keying (CSK) and the recei vers count the number of molecules at the time of sampling. As customary , we set the sampling time so that the number of molecules at receiv ers for the corresponding transmitters is maximized. As sho wn in Fig. 2, the channel has memory and due to the inter -symbol- interference (ISI), the recei ver counts the molecules from previous samples of the corresponding transmitter . Similarly , the molecules from the current and previous samples of the non-corresponding transmitters are known as inter-link- interference (ILI). Number of channel taps (L) for each link is related to the system geometry , configuration and sample interval. Specifically , we can eliminate the ISI and ILI by making the sample interval large enough and putting each pair of transceivers far enough from the other pairs. Howe ver , this case is not considered due to the demands for high data rate per unit of space. Hence, we ha ve to face the ISI and ILI in the MC system and try to mitigate their effects. The observed number of molecules at sampling time k and recei ver j is y j [ k ] = M X i =1 L − 1 X ` =0 c ij [ `, k ] x i [ k − ` ] + v j [ k ] (1) where L is the memory taps, c ij [ `, k ] is a random variable and denotes to the number of molecules observed at the receiv er j from transmitter i due to the release of N molecules at the time interval [ k − ` ] . Case i = j refers to the paired transmitter- receiv er, otherwise it refers to the inter-link interference. x i [ k ] ∈ { 0 , 1 } is the transmitted symbol at the time interval k from transmitter i . The number of molecules c ij [ `, k ] can be approximated as a Poisson random variable with a mean value ¯ c ij [ ` ] : c ij [ `, k ] ∼ P oiss ( ¯ c ij [ ` ]) . Additionally , v j [ k ] is the number of external noise molecules detected at the recei ver j at time interval k . Noise molecules could originate from the remaining channel taps from all transmitters not considered in model, and any external source. Hence, we can consider the noise as a Poisson with a mean ¯ v j : v j [ k ] ∼ P oiss ( ¯ v j ) [16]. Fig. 2. Impulse response, ¯ c ij ( t ) , of a 2 × 2 D-MIMO system at Rx 1 vs time, for 3 emissions of molecules with time spacing 0 . 2 ms : ILI and ISI are black dots. I I I . P R O B L E M D E FI N I T I O N Assume that x i = [ x i [1] , x i [2] , ..., x i [ K ]] T is a binary training sequence with length K for transmitter i . T o avoid edge effect due to the ISI, in CIR estimation we employ y j [ k ] for k ≥ L . Therefore, the K − L + 1 samples are used for CIR estimation of the j -th receiv er . The mean number of molecules at receiv er j at time k is ¯ y j [ k ] = E { y j [ k ] } = M X i =1 L − 1 X ` =0 ¯ c ij [ ` ] x i [ k − ` ] + ¯ v j (2) which L ≤ k ≤ K . Eqn. (2) can be written compactly as: ¯ y j [ k ] = X T [ k ] ¯ C j (3) where the following notations are used: x [ k ] = [ x 1 [ k ] , x 2 [ k ] , ..., x M [ k ]] T X [ k ] = [ x T [ k ] , x T [ k − 1] , ...., x T [ k − L + 1] , 1] T ¯ c j [ ` ] = [¯ c 1 j [ ` ] , ¯ c 2 j [ ` ] , ..., ¯ c M j [ ` ]] T ¯ C j = [ ¯ c T j [0] , ¯ c T j [1] , ..., ¯ c T j [ L − 1] , ¯ v j ] T here, x [ k ] is a M × 1 vector denotes to the training sequence of all transmitters at time k , ¯ c j [ ` ] is a vector with dimension M × 1 denoting the expected number molecules at receiver j at time k due to the transmission of N molecules at time k − ` from transmitter T x i , and x [ k − ` ] is its corresponding training sequence. Additionally , ¯ C j is a ( M L + 1) × 1 vector gathering all channel memory taps of receiv er j due to the ISI and ILI from all transmitters and also noise ¯ v j , and X [ k ] with dimension ( M L + 1) × 1 is its corresponding training sequence vector at time interv al k . The expected number of molecules at receiv er j during the time interv als L ≤ k ≤ K is defined as ¯ y j = [ ¯ y j [ L ] , ¯ y j [ L + 1] ..., ¯ y j [ K ]] T , and the D-MIMO relation for receiver j is written ¯ y j ( K − L +1) × 1 = X T ( K − L +1) × ( M L +1) ¯ C j ( M L +1) × 1 (4) where X is a ( M L + 1) × ( K − L + 1) con volution matrix of training sequences including memory of pre vious samples due to the channel taps, and it is defined as X = [ X [ L ] , X [ L + 1] , ..., X [ K ]] (5) Finally , we define ¯ Y = [ ¯ y 1 , ¯ y 2 ..., ¯ y M ] and we compactly write the global D-MIMO relation into ¯ Y ( K − L +1) × ( M ) = X T ( K − L +1) × ( M L +1) ¯ C ( M L +1) × M (6) where ¯ C is the ( M L + 1) × M global channel matrix, and it is defined as ¯ C = [ ¯ C 1 , ¯ C 2 , ..., ¯ C M ] . (7) The matrix of all the observed number of molecules at the M receiv ers contain the Poisson random variables with mean equal to ¯ Y : Y = P oiss ( ¯ Y ) (8) which means each entry of the observed matrix Y , is Poisson random v ariable with mean equal to the corresponding entry of ¯ Y . The probability density function (PDF) of all observations at all recei vers are the product of the Poisson distribution of each observation at each receiver f Y ( Y | ¯ C , X ) = K Y k = L M Y j =1  X T [ k ] ¯ C j ) y j [ k ] exp( − X T [ k ] ¯ C j ) y j [ k ] ! (9) According to (4), we can analyze the performance of each receiv er independently to make sure that all M receiv ers are simultaneously working optimally . Therefore, the PDF of the observations of j -th receiv er is f y j ( y j | ¯ C j , X ) = K Y k = L  X T [ k ] ¯ C j ) y j [ k ] exp( − X T [ k ] ¯ C j ) y j [ k ] ! , (10) for maximum likelihood estimation belo w . I V . D - M I M O C H A N N E L E S T I M AT I O N In this section sub-optimal maximum lik elihood (ML) and least squares (LS) estimators for the D-MIMO system are deriv ed. Then, CRB for each recei ver is computed. A. Maximum Likelihood estimator Maximum likelihood (ML) D-MIMO CIR estimator finds the CIR which maximize the likelihood of observation vector y j ˆ ¯ C ML j = argmax ¯ C j ≥ 0 f y j ( y j | ¯ C j , X ) = argmax ¯ C j ≥ 0 L y j ( y j | ¯ C j , X ) (11) where the log likelihood function and is gi ven by L y j ( y j | ¯ C j , X ) = K X k = L h − X T [ k ] ¯ C j + y j [ k ] ln( X T [ k ] ¯ C j ) i (12) The ML estimate of the CIR for the D-MIMO channel at receiver j is obtained by solving a system of non-linear equations giv en below [14]: K X k = L h y j [ k ] X [ k ] X T [ k ] ¯ C j − X [ k ] i = 0 (13) Remar k 1 : W e note entries of ¯ C j are positive semidefinite. Howe ver , ML estimator could estimate a negati ve value for some elements of the ¯ C j . Sub-optimal solution is to set to zero all the negativ e entries of the estimated CIR. This heuristic approach was adopted for single link MC [14], and showed therein a negligible loss of performances compared to the optimal ML. Therefore, sub-optimal solution of (13) is highly preferred in D-MIMO channels due to its simplicity . B. Least Squar es Estimator The least squares (LS) method chooses ¯ C which minimizes the sum of the square errors at all recei ver from the observ ation vector Y , ˆ ¯ C LS = argmin ¯ C ≥ 0 k  k 2 . (14) where  = Y − E { Y } = Y − X T ¯ C . The LS estimate of the CIR for D-MIMO channel is ˆ ¯ C LS = h ( X X T ) − 1 X Y i . (15) Minimization of (14) is a constrained optimization problem with constraint C ≥ 0 for entries. In case there exist a stationary point, this is the global optimum solution. In case the stationary point does not exist, sub-optimal solution is to set all negati ve elements of C to zero. Optimal solution for (14) is introduced in [14], and the authors showed that for K lar ge, there e xist a stationary point, and for small lengths, the performance loss is very negligible. Again, we prefer the sub-optimal solution for D-MIMO system due to its simplicity . C. Cramér -Rao Bound The Cramér-Rao bound (CRB) sets the lo wer bound on the cov ariance of any unbiased estimator of a deterministic parameters. Let ˆ ¯ C j be the unbiased estimator of ¯ C j , the CRB sets the bound of the cov ariance cov ( ˆ ¯ C j )  I − 1 ( ¯ C j ) (16) where I ( ¯ C j ) is the Fisher information matrix of ¯ C j and is giv en by I ( ¯ C j ) = E y j {−∇ 2 ¯ C j ¯ C j L y j ( y j | ¯ C j , X ) } (17) Therefore, the CRB at receiv er j is given by C R B j = tr { I − 1 ( ¯ C j ) } = tr (" K X k = L X [ k ] X T [ k ] X T [ k ] ¯ C j # − 1 ) . (18) Notice that (16) implies that the difference cov ( ˆ ¯ C j ) − I − 1 ( ¯ C j )  0 is positi ve semidefinite. Fig. 3. T opological model for a 2 × 2 D-MIMO system. Both transmitters use the same information molecules (pentagon shapes), and molecules here hav e different colors according to the corresponding transmitter . V . T R A I N I N G S E Q U E N C E D E S I G N In this section, we present a method for designing the training sequences for estimating the CIR of a D-MIMO channel. As shown in (18), the CRB for a giv en system is a function of training sequences. Therefore, we can find a set of training sequences that minimize the CRB of all receiv ers. In other words, the CRB of a specific receiver depends on the training sequence of all transmitters which hav e interference with it. In general, for a M × M D-MIMO system, we have to design M different training sequences to minimize the CRB of all receivers simultaneously . Ho wev er , in practice we do not need to design M training sequences, because ILI for f ar transmitters is negligible and thus we neglect their interference channels but consider them as an augmented noise source in v j . In order to find a suitable set of training sequences that simultaneously minimize all CRBs, we consider following constraints: 1) the training sequences should be molecularly efficient by minimizing the fraction of molecules used for channel estimation, and 2) transmitters can not be silent for many consequent interv als. In detail, for a training sequence of length K, we consider sequences with maximum K / 2 ones, consequently transmitting maximum N K / 2 molecules, and the maximum consequent zeros are considered 4 time intervals. X is the sets of all possible training sequences that meet the above criteria. [ x 1 , x 2 , ..., x M ] = argmin x i ∈X { C R B 1 , . . . C RB M } (19) Remark 2: Accuracy of CIR estimation depends on the training sequence length, hence K should chosen carefully . For large K , it is difficult to search among all suitable sets and find the optimum ones. Therefore, we look for an optimum training sequence with smaller length K 1 , and we concatenate it to b uild a longer training sequence of length K . W e observed that if K 1 is wisely selected, concatenating w ould not impair the performances. Fig. 4. Comparison of ML and LS estimators to CRB in terms of MSE in dB vs. the training sequence length K with L = 3 . V I . P E R F O R M A N C E E V A L U A T I O N O F T H E D E S I G N E D D - M I M O S Y S T E M In this section, we present a 2 × 2 D-MIMO configuration and we compare the performances of the estimators introduced in this paper for ON-OFF keying signaling. W e have generated the CIR according to the analytical models proposed in [8], [5] for MC systems. Ho wev er , there is no constraint in the value of CIR to be estimated. Dif fusion coefficient value is 10 − 9 m 2 /s and it is compatible to the normal values of diffusion of most of molecules in water at room temperature. Choosing bit interval time is a trade-of f between bit rate and total number LM of ISI and ILI of the channel memory to be estimated. Bit interval time is T int = 0 . 2 ms , all transmitters release N = 10 5 molecules and the receiv ers counts once the number of molecules per each symbol at time which CIR of the pair transmitter is expected to be maximum. For simplicity , the mean of noise is chosen as ¯ v j = 0 . 3 ¯ c j j (0) . The number of channel taps for both ISI and ILI link are considered L = 3 , so ¯ c ij [ L ] ≤ 0 . 05 ¯ c j j [0] . The 2 × 2 D-MIMO system is shown in Fig. 3. W e hav e assumed that the distance between transmitter and receiv er is d = 400 nm and the inter-distance is h = 200 nm . Spherical receiv er with radius 50 nm is assumed. Positions of the mentioned entities are fluctuating: P T x 1 = (0 + δ x 1 , 0 + δ y 1 , 0 + δ z 1 ) , P T x 2 = (0 + δ x 2 , h + δ y 2 , 0 + δ z 2 ) , P Rx 1 = ( d + δ x 3 , 0 + δy 3 , 0 + δz 3 ) , P Rx 2 = ( d + δ x 4 , h + δy 4 , 0 + δz 4 ) , with δ x,y ,z ∼ N (0 , σ 2 ) , and σ 2 = 50 nm . Since the T x and Rx are not fixed in the position, the channel is varying in time. While the entities are fixed, the CIR at the receiv ers are: ¯ C 1 = [60 . 21 , 41 . 58 , 9 . 11 , 8 . 71 , 3 . 83 , 3 . 74 , 18 . 06] T and ¯ C 2 = [41 . 58 , 60 . 21 , 8 . 71 , 9 . 11 , 3 . 74 , 3 . 83 , 18 . 06] T . W e note that in each realization the CIR is dif ferent, ¯ C 1 6 = ¯ C 2 , because the position of transmitters and receiv ers are changing with normal distrib ution with σ 2 = 50 nm . The training sequences are designed according to (19) with the length K 1 = 16 and they are x 1 = [1 , 1 , 1 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0 , 1 , 1 , 0 , 0 , 1] T and x 2 = [1 , 1 , 1 , 0 , 1 , 0 , 0 , 0 , 1 , 1 , 1 , 0 , 0 , 0 , 0 , 1] T . Longer training sequences are constructed by concatenating these training sequences as detailed in section 5. In this problem, each receiv er has to estimate LM + 1 = 7 variables, for a total of 14 variables. The results in Fig. 4, are Monte Carlo simulations with 1000 random CIR realizations. In terms of mean square error (MSE), E  ˆ ¯ c j − ¯ c j || 2  in dB vs. the training sequence length (K) for the ML and LS (Fig. 4) estimators, respecti vely . The MSE decreases with increasing the training sequence length as we expected. W e can notice that training sequences are designed such that both recei vers hav e optimum performances for both estimators as they attain the corresponding CRB. Fig. 5, shows the Normalized MSE which is defined by M S E N j = E  || ˆ ¯ c j − ¯ c j || 2  || E { ¯ c j }|| 2 . (20) The value of the normalized MSE is much lower , around 38dB, than the MSE. As we can see in Fig. 5, the perfor- mance of the ML estimator outperforms the LS estimator by approximately 1 dB. Howe ver , the LS estimator is preferred due to its simplicity respect to the ML estimator, because our bio-based receiv ers hav e limited computational capabilities. In applications where receivers send the data to the external computers, the ML estimator is preferred because it reaches to the CRB bound. V I I . C O N C L U S I O N S In this paper , we hav e presented a training-based channel es- timation for M × M D-MIMO system. W e hav e dev eloped ML and LS estimators which consider ISI of the pair transmitter - receiv er and also inter-link diffusi ve interference. W e have also derived the Cramér-Rao bound for each receiv er . T raining sequences have been designed such that minimize the CRB for all receiv ers simultaneously . The ML estimator outperforms the LS estimator and the corresponding MSE reaches the CRB. Howe ver , ML estimator needs to solve a system of non-linear equations and this makes it computationally not preferable. In the other hand, LS estimator is very simple to design at the price of a small performance degradation (1 dB for the example considered here) compared to the CRB. R E F E R E N C E S [1] T . Nakano, T . Suda, M. Moore, R. Egashira, A. Enomoto, and K. Arima, “Molecular communication for nanomachines using intercellular calcium signaling, ” in Nanotechnology , 2005. 5th IEEE Conference on , pp. 478– 481, IEEE, 2005. [2] N. Farsad, H. B. Y ilmaz, A. Eckford, C.-B. Chae, and W . Guo, “ A com- prehensiv e survey of recent adv ancements in molecular communication, ” IEEE Communications Surveys & T utorials , v ol. 18, no. 3, pp. 1887– 1919, 2016. [3] I. Akyildiz, M. Pierobon, S. Balasubramaniam, and Y . K oucheryavy , “The internet of bio-nano things, ” IEEE Communications Magazine , vol. 53, no. 3, pp. 32–40, 2015. [4] T . Nakano, T . Suda, Y . Okaie, M. J. Moore, and A. V . V asilakos, “Molecular communication among biological nanomachines: A layered architecture and research issues, ” IEEE transactions on nanobioscience , vol. 13, no. 3, pp. 169–197, 2014. [5] S. Kadloor, R. S. Adve, and A. W . Eckford, “Molecular communication using brownian motion with drift, ” IEEE T ransactions on NanoBio- science , vol. 11, no. 2, pp. 89–99, 2012. Fig. 5. Normalized MSE in dB vs. training sequence (K) for ML and LS estimators with L = 3 [6] T . Nakano, A. W . Eckford, and T . Haraguchi, Molecular communication . Cambridge Univ ersity Press, 2013. [7] M. S. Kuran, H. B. Y ilmaz, T . T ugcu, and I. F . Akyildiz, “Modulation techniques for communication via diffusion in nanonetworks, ” in Com- munications (ICC), 2011 IEEE International Confer ence on , pp. 1–5, IEEE, 2011. [8] A. Noel, K. C. Cheung, and R. Schober, “Optimal receiver design for diffusi ve molecular communication with flow and additive noise, ” IEEE transactions on nanobioscience , vol. 13, no. 3, pp. 350–362, 2014. [9] L. C. Cobo and I. F . Akyildiz, “Bacteria-based communication in nanonetworks, ” Nano Communication Networks , vol. 1, no. 4, pp. 244– 256, 2010. [10] N. Farsad, Y . Murin, A. Eckford, and A. Goldsmith, “On the capacity of diffusion-based molecular timing channels, ” in Information Theory (ISIT), 2016 IEEE International Symposium on , pp. 1023–1027, IEEE, 2016. [11] N.-R. Kim, A. W . Eckford, and C.-B. Chae, “Symbol interval opti- mization for molecular communication with drift, ” IEEE tr ansactions on nanobioscience , vol. 13, no. 3, pp. 223–229, 2014. [12] B.-H. Koo, C. Lee, H. B. Yilmaz, N. F arsad, A. Eckford, and C.-B. Chae, “Molecular mimo: From theory to prototype, ” IEEE J ournal on Selected Ar eas in Communications , vol. 34, no. 3, pp. 600–614, 2016. [13] L.-S. Meng, P .-C. Y eh, K.-C. Chen, and I. F . Akyildiz, “Mimo com- munications based on molecular diffusion, ” in Global Communications Confer ence (GLOBECOM), 2012 IEEE , pp. 5380–5385, IEEE, 2012. [14] V . Jamali, A. Ahmadzadeh, C. Jardin, H. Sticht, and R. Schober, “Channel estimation for diffusi ve molecular communications, ” IEEE T ransactions on Communications , vol. 64, no. 10, pp. 4238–4252, 2016. [15] H. Arjmandi, A. Gohari, M. N. Kenari, and F . Bateni, “Diffusion- based nanonetworking: A new modulation technique and performance analysis, ” IEEE Communications Letters , vol. 17, no. 4, pp. 645–648, 2013. [16] A. Noel, K. C. Cheung, and R. Schober , “ A unifying model for external noise sources and isi in diffusi ve molecular communication, ” IEEE Journal on Selected Ar eas in Communications , vol. 32, no. 12, pp. 2330– 2343, 2014.