GSM based CommSense system to measure and estimate environmental changes
GGSM based CommSense system to measure and estimateenvironmental changes
Abhishek Bhatta
Electrical Engineering DepartmentUniversity of Cape TownEmail: [email protected]
Amit Kumar Mishra
Electrical Engineering DepartmentUniversity of Cape TownEmail: [email protected] 15, 2018
Abstract
Facilitating the coexistence of radar systems with communication systems has been a major area ofresearch in radar engineering. The current work presents a new way to sense the environment using thechannel equalization block of existing communication systems. We have named this system CommSense.In the current paper we demonstrate the feasibility of the system using Global System for Mobile Communications(GSM) signals. The implementation has been done using open-source Software Defined Radio (SDR)environment. In the preliminary results obtained in our work we show that it is possible to distinguishenvironmental changes using the proposed system. The major advantage of the system is that it is inexpensiveas channel estimation is an inherent block in any communication system and hence the added cost to makeit work as an environment sensor is minimal. The major challenge, on which we are continuing our work, ishow to characterize the features in the environmental changes. This is an acute challenge given the fact thatthe bandwidth available is narrow and the system is inherently a forward looking radar. However the initialresults, as shown in this paper, are encouraging and we intend to use an application specific instrumentation(ASIN) scheme to distinguish the environmental changes.
Classical radar systems have been designed primarily for military operations. Of late many interesting ways ofusing radio-frequency spectrum for radar purpose have been under research. One such concept is commensal The word commensal has been borrowed from biology in which this represents co-existence of two species out of which one isbenefited and the other remains unaffected. a r X i v : . [ c s . OH ] M a y r passive radar, which uses signals of opportunity to detect targets without affecting the functionality of theparent system [6, 12–14, 26, 28, 29].One implementation of commensal system is the recently active area of Passive Bistatic Radar (PBR)[7, 9, 11, 16, 23].In our work we shall describe a kind of commensal radar which we call communication sensor (CommSense)[17, 18] system. It uses the channel estimation processing in communication systems to estimate changes inthe environment. This system can, potentially, be used to monitor land terrain, sea state or even naturaldisasters. The novelties of the proposed system are as follows.First of all, this system is built upon existing communication systems. The cost of implementing such asystem is therefore considerably low. Secondly, unlike PBR systems it does not process the informationusing correlation, rather it uses a known training sequence to estimate changes in the environment. Ittherefore eradicates the need for a reference antenna pointing towards the transmitter. We have implementedthe prototype system using open source GNURadio software and Nuand BladeRF × SDR hardware [8,20].It must, however, be noted that there are two major challenges in implementing this system. Firstof all, the conventional concepts of resolution will not be valid here mainly because of the fact that thebandwidth available to us is very narrow. This makes the measurement system a severely ill-posed inverseproblem. Secondly, this is a forward looking and non-coherent radar system which theoretically limits theamount of information it can capture. We hypothesize that these two problems can be solved using theapplication specific instrumentation framework [19, 24]. This will be our future work, whereas in this paperwe demonstrate the feasibility of the concept of sensing the environment from communication signals.The rest of the paper is organized as follows. Section 2 gives a basic understanding of GSM protocol andan introduction to the design of this system. Section 3 provides details about the real time implementationof the system. Section 4 shows how the data analysis is done, specifically focussing on the ProbabilityDistribution Function (PDF) analysis and clustering. Finally Section 5 concludes the work and shows thescope for future developments.
Any signal transmitted over the channel gets affected by the channel itself. This effect is minimized inwireless communication systems using post-reception processing blocks called channel estimation and equalization[21, 22]. The hypothesis of our work is that the channel information from the received signal can beharnessed and used to statistically monitor the changes in the environment. For example we will be ableto differentiate between different weather conditions or produce an alternate model (which, potentially, canTHIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272ive more information) for sea waves and mountainous terrains. We call this system CommSense system.The reader is suggested to refer [18] for system level information of a generic CommSense system.The communication technology infrastructure used for the current implementation of the CommSense isGSM. This is because it is currently the most widely used wireless communication technology in SouthAfrica. We are also working towards investigating its feasibility using LTE infrastructure [25]. GSMprovides coverage of almost the entire nation. It is a wireless communication protocol developed by EuropeanTelecommunications Standard Institute (ETSI) to provide a standard in the wireless telecommunication industry.CommSense system relies on the broadcast signals transmitted by the base stations for phones to detect andconnect to the particular base station, thus not hampering the communication system in any way. GSMtransmits a known training sequence in every frame so that the receiver is able to detect and cancel thechannel interference. Using the information from this, the channel states can be estimated.Figure 1 shows the concept diagram of the application. Signal from the base station reaches the mobilestation through multiple different paths and is affected by the physical properties of each path. Out of all thesereturns, the signal with highest Signal to Noise Ratio (SNR) is extracted by the method of preliminary channelequalization and used for the purpose of communication. In our work we shall be using this preliminarychannel equalization block. It, however, can be noted that there are further blocks in the processing chainin the receiver which use a range of algorithms to take care of channel noise. For the purpose of our workwe shall use channel equalization to mean this preliminary channel equalization processing. In order toimplement channel equalization the communication system transmits a known sequence of bits with everyframe. These known bits, also known as training sequence, are extracted from the channel convolutedreceived signal and used to extract The difference between the known and expected signal. This differencebetween the known and the received signal gives the state of the channel through which the signal hastravelled. In the proposed system this estimated information is used to find the changes in the channel statein different environmental conditions.
CommSense system acts as a commensal radar, eliminating the need of a new transmitter. As the processingpower of smart-phones has significantly improved, it is planned to eventually implement CommSense as amobile application. For the purposes of proving the hypothesis it has been implemented on the SDR platformwith open source software. The hardware used for the implementation is BladeRF × manufactured byNuand [20]. It is a Universal Serial Bus (USB) powered SDR board with an on-board transmitter and receiverport. It can operate over a frequency range of 300 MHz to 3.8 GHz which gives the opportunity to implementTHIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272 irect Path R e fl e c t e d p a t h R e fl e c t e d p a t h BTS
Mobile Station
Channel Estimation
Real time plot of the channel e ff ects Sea Waves
Ground Surface R e fl e c t e d p a t h R e fl e c t e d p a t h Clouds
Figure 1: System Overview.the proposed system using GSM technology. The software used for implementing the real time channelestimation system is GNU Radio [8]. It is an open source software tool-kit, distributed under GNU generalpublic license, that provides the opportunity to implement signal processing blocks using open source RFhardware for research purposes.GSM is a wireless communication protocol that acts as a host to this system [1, 2, 10]. It uses FrequencyDivision Multiple Access (FDMA) and Time Division Multiple Access (TDMA) in order to accommodatemultiple users. The various frequency bands of operation for this technology are 800 MHz, 900 MHz, 1800MHz, 1900 MHz. GSM works in frequency duplex mode thus having a different frequency for uplink anddownlink. GSM900, used here for data collection, uses 880-915 MHz for uplink and 925-960 MHz fordownlink communications. These bands are further divided into 200 kHz bands and are separated by anumber called Absolute Radio Frequency Channel Number (ARFCN). Each wireless channel is divided intoTDMA frames with 8 time slots of 577 µ s each. Each timeslot is time shared between the mobile and thebase station. The logical channels are piggybacked on the physical channel, described in [1].Each frame in the GSM system transmits a known training sequence for the receiver to detect and reducethe possibility of transmission error. This functionality provided by GSM is used in this project to estimatethe channel effects on the signal and analyse it to find the differences in the various types of terrain andenvironmental conditions. A frame in GSM system is shown in Figure 3. This structure shows the normalburst that has been used for the extraction of the parameters. This burst has 114 bits of information transmittedwith tail bits, flags and a guard period to identify the start and stop of the burst and different datasets. Thetraining sequence is 26 bits long out of which the central 16 bits are used for channel equalization. In theTHIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272urrent implementation these equalized channel values are used to determine the surrounding environment.A guard period is added at the end of each frame to provide a window of error against distortions that occurdue to the rise and fall time of the signal. Figure 2 shows the block diagram of the implementation. The received GSM signal can be mathematicallyrepresented as equation (1). Here, the received signal y ( t ) is the convolution of the transmitted signal x ( t ) and Channel Impulse Response (CIR) h ( t ) in the presence of Additive white Gaussian noise (AWGN) n ( t ) .The transmitter sends a known training sequence in each frame as shown in Figure 3, which is divided intoreference length of P and guard period of L bits [22]. This equation can be represented in the matrix formand is shown in equation (2). y ( t ) = h ( t ) ∗ x ( t ) + n ( t ) (1) y = Mh + n (2)In equation (2) M represents a Toeplitz like matrix structure of the given training sequence as shownin equation (4). This matrix is made of the m array that is shown in equation (3), which is actually theoversampled information from each received frame after filtering out the Gaussian noise. m = [ m m .... m P + L − ] T (3) M P + Cl − × Cl = m · · · m m · · · m m m · · · ... ... ... ... ... m P m P − m P − · · · m m P m P − · · · m m P · · · m ... ... ... ... ... m P (4) CIR is a black-box linear model of the channel effects.
THIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272 ignal Source x(t)
Channel + Noise h n(t)
Receiver FilterChannel Estimator MLSE Detector a y(t) y(t) h est est Communication Applications
Figure 2: Channel Estimation Block Diagram. h = [ h h h · · · h n ] T (5)In equation (4), the matrix M is a P + Cl − × Cl matrix where P is the reference length of trainingsequence and Cl is the length of the CIR. equation (5) shows the CIR array, here each value of the arrayrepresents the channel parameters of one reflected signal. Thus the length of the array determines how manymulti-path signal are being analysed. The Least Squares algorithm finds the CIR by minimising the squarederror quantity which in the presence of white Gaussian noise gives equation (6). h est = ( M H M ) − M H y (6)Here, M H denotes Hermitian transpose matrix and () − denotes matrix inverse. In equation (6) we canobserve that the received signal matrix y is multiplied to a matrix also known as the pseudo inverse matrix of M . Figure 4 shows the plot for simulated CIR in the presence of AWGN. This result is for a single GSMframe being transmitted over a known channel of length 5 taps. The estimated value is very similar to theoriginal CIR value which shows that the algorithm is working properly. Then this is implemented in real-time Burst
Tail bits Information Flag FlagTraining Sequence Guard PeriodInformation Tail bits
Figure 3: GSM Normal Channel Frame Structure.THIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272s shown in Figure 5.
Filter Taps A m p li t u d e Channel ImpairementsEstimated Channel
Figure 4: Simulated Least Square Channel Estimation in presence of AWGN.
The system is implemented in GNU Radio, which does the channel estimation and plots the CIR in real timeas shown in Figure 6. The hardware used is BladeRF × , which works at a frequency band of 300 MHz to3.8 GHz with a maximum physical bandwidth of 28 MHz. The antenna used for this specific implementationis an off the shelf quad band GSM antenna. The accuracy of the system is proved by decoding the receivedbit-stream and analysing it in Wireshark (a packet analyser tool), which gives the basic information as thelocation of the base station, carrier service provider and so on. The steps taken to implement the channel equalizer in real-time are mentioned here. Correlation between thereceived signal and the known training sequence is used to extract the differences between the two signalsand this information is saved in a file for further analysis.The SDR hardware, bladeRF, receives the GSM signal downconverts it according to its Phase LockedLoop (PLL) clock frequency and gives out the In-phase ( I ) and the Quadrature ( Q ) components of thesignal. The receiver first needs to get synchronized with the base station, the Synchronization Channel (SCH)burst is used for this purpose. The SCH burst consists of a 64 bit long training sequence and transmits theTHIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272ame sequence in every timeslot allocated to it, thus it is easily distinguishable and the receiver can getsynchronized. Once the receiver synchronizes with the GSM base station, the frequency offset is calculatedusing the Frequency Correction channel (FCCH).The received burst type is identified in two steps first with the help of the carrier index of the burst andthen with the burst number. Carrier index zero (C0) is a special case where only Broadcast Channel (BCCH)data is transmitted. BCCH carries a repeating pattern of system information such as identity, configurationand available features for the base station. In case of other carrier indices the burst type is determined by theburst number. The most commonly transmitted burst is the normal burst thus in this implementation we arefocussing on extracting the channel information from the normal burst.In order to receive a burst completely the receiver needs to wait for two consecutive guard periods.Although there will be an overlap between the guard periods of two consecutive bursts, this is necessaryto ensure full burst is received. Once a normal burst is identified the start and the stop position of the burstis determined eliminating the guard period as it consists of 8.25 bits on each side and is distinguishable.The central 26 bits of the remaining bits are extracted as it constitutes of the received training sequence.The training sequence extracted here is a complex data so the already available sequence also needs to beconverted into a complex number. The two complex training sequences are correlated with each other andthe correlation information is saved in a buffer. A pointer is placed at the beginning of the buffer and basedon predefined CIR length the pointer is moved from the beginning to the length of the CIR. The valuescorresponding from the beginning to the CIR length is saved in an array and then moved to a file where itis saved for further analysis. In case of communication systems the maximum absolute value of the CIR ischosen discarding the rest of the values and that information is used to equalize the channel and get bettercommunication. For the purpose of this implementation the entire channel information is necessary thus allof it is saved.In order to check if the receiver is working as intended the entire receive structure is made referring tothe gr-gsm libraries available for GNU Radio and the bursts are printed out using a message printer. Figure 6shows the plot of the entire implementation including the GSM frequency spectrum, received bursts and thereal time estimated channel values. Figure 2 shows the entire flow of the implementation, where the signal source is the transmitting base stationand the processing is done on over the air signals. The receiver filter block removes most of the noise fromthe signal and shows the GSM frequency spectrum as shown in Figure 6. The data is used to estimate theTHIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272
100 200 300 400 500 600
Channel Characterization Values A m p li t u d e Time Domain Absolute Channel Impulse Response
Figure 5: Estimated Normal channel absolute values in terms of the I and Q component, in this figure isrepresented by h est = (cid:112) I + Q . Real Time Estimated Channel ValuesGSM Frequency SpectrumReceived bursts BladeRF SDR Hardware
Figure 6: GNU Radio implementation of Channel Equaliaztion with BladeRF × .channel and store the values in a file for further analysis. The darker blocks in Figure 2 are the normal GSMprocessing blocks which are not affected by this system as it runs parallel to the communication system.One of the major challenges faced during the implementation is finding the start-point of the burst. Thisis done with the help of the guard period. The receiver waits for two consecutive guard periods, once receivedit considers the start point of the burst. In the case of an error in receiving the guard periods, the system dropsthe current set of data and waits for the next set and then repeats until it gets two consecutive guard periods,thereby eliminating the chances of any ambiguous data.THIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272 Analysis of Captured Data
To prove the hypothesis the real-time channel information is captured in different locations and classified asStatistical Characterization, Distinguishing Events and Special case. Table 1 lists details of all the datasetscollected and analysed in this work. The estimated channel information is analysed to show the changesin the channel state in different environmental condition. Three different types of analysis are presentedin this section, two in statistical domain and one in data/time domain. In the first approach we estimatethe PDF of the channel values captured at different scenarios. This shows that the estimated PDF’s varyas the environmental conditions change which means that a “hypothesis test” kind algorithm can be usedto distinguish different environment types from the estimated channel values. In rest of the paper we shalluse simply PDF to mean estimated PDF. In the second approach we perform Chi-square test to assess thegoodness of fit between the observed data and the theoretically expected dataset. In the third approachwe check for the clustering of the estimated channel values for different capture scenarios. For easiervisualisation we used principal component analysis (PCA) to reduce the dimension of the data. This showsthat the data are highly clustered in PC domain.
Captured DataType Location CommentsStatistical Characterization
Sea Beach To check the effects of the sea waves on the signal,multiple locations may also be testedHill next to a road Check the effects on the signal in the presence of ahill and see repeatability (safe location still needs to bedecided)Highland Captures taken at the road next to a steep terrainHeavy Rain, slightlyhumid climate, hot day Captures taken at the same location at different instancesof the day showing the differences due to differentweather conditionsParking space Captures are taken on different times of the day. Oneset showing when the place is full and the other showingempty parking lot
Distinguishing events
Train Station At a train station with and without a train in proximityJ stairs At the stairs to the entrance of an approximately 10 metertall building.Bus With and without the presence of a bus in the proximity,two different sets of captures were taken at differentlocations
Special case
Car Captures taken with and without the presence of a car inproximity of the receiverCorner Reflector Captures are taken when a corner reflector is placed at adistance of 2 meter from the antenna in various differentconfigurations ’V’ is vertical, ’H’ is horizontal, ’VH’refers to a dihedral corner reflector
Table 1: The different sets of captured data.THIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272 .1 PDF Analysis
In this set of analysis we plot the empirical PDF of the data collected from different scenarios as describedin Table 1. In addition to the empirical PDF we also plot the PDF of the data for four different theoreticalPDF models, i.e. Rayleigh, Normal, Log-normal and Gamma distributions. A brief description of thesedistributions and maximum likelihood estimation (MLE) expression of their parameters are given in AppendixI. Figure 7a and 7b show the distribution plot of two different datasets in similar conditions. Both thecaptures are taken near the entrance of a building about 10 meters tall. The difference between the twocaptures is the wind speed. During the capture of Figure 7a the wind speed was slower compared to thecapture of Figure 7b. The location picture of the J stairs capture set is provided in Appendix II.The plots display index 12 of two different datasets whose indices go up to a total of 40 which is thenumber of multipath signals extracted from the received signal. Here we observe the similarities in thedistribution pattern of the captures. Although the empirical PDF are not completely matched, the generalpattern is similar.This proves that the data captured from a single location will generally follow a similar pattern in aparticular environmental condition (wind, humidity, temperature, etc.), thus aiding in providing enoughinformation to map the surface and the environment at the time of the capture.Next we analyse the empirical PDF and theoretical PDF of estimated channel values from differentscenarios. Some of the example empirical PDFs are shown in Figures 8,9,10,11. Each curve shows theempirical PDF of a particular case. Keeping the y -axis constant shows the differences in the given datasetsclearly. The datasets presented here include the case of car in proximity, near a big building, bus in proximity,full parking space, train in proximity and the special case with and without the corner reflector. Dataset Shape Mean Variance Skew KurtosisLog-normal Gammajameson stairs .
762 1 .
61 1 . e − . e − .
15 6 . jameson stairs2 .
767 1 .
58 1 . e − . e − .
38 7 . rts t .
28 4 .
80 4 . e − . e − .
58 0 . roof bcr loc2 .
74 1 .
67 2 . e − . e − .
95 4 . roof vh metal cr .
46 2 .
65 1 . e − . e − .
16 1 . north stop jammie loc1 .
63 2 .
15 1 . e − . e − .
95 5 . car .
40 2 .
91 1 . e − . e − .
97 1 . rhodes mem .
29 5 .
55 9 . e − . e − .
95 1 . parking full .
43 3 .
07 1 . e − . e − .
14 1 . Table 2: Moments of captured datasets.In Figure 8 and 9 we observe the change in variance with change in the target conditions. In the caseof a bus, the distribution has a lower variance which increases as the target changes to a car and even moreTHIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272 .00 0.02 0.04 0.06 0.08 0.10 0.12
Estimated Channel Values -1 P r o b a b ili t y D i s t r i b u t i o n Distribution Fitting J_Stairs rayleighnormlognormgamma (a) Empirical PDF near a 10 meter tall building in the presence ofslow wind speed.
Estimated Channel Values -1 P r o b a b ili t y D i s t r i b u t i o n Distribution Fitting J_Stairs2 rayleighnormlognormgamma (b) Empirical PDF near a 10 meter tall building when the windspeed is high.
Figure 7: Empirical PDF comparison of th index taken from a 40 index dataset near a building withdifferent wind speeds. In the legend, rayleigh , norm , lognorm and gamma represents theoreticalPDF assuming Rayleigh distribution, Gaussian distribution, Lognormal distribution and Gamma distributionrespectively.in case of the train. Looking at Figure 8b and 9b a similar pattern can be observed with minimum changesin the distribution. This is because both the captures represent cars in a flat space, the only difference is inFigure 8b there is only one car in proximity and in Figure 9b there are many cars which creates the observabledifferences in the plots. Figure 10 gives the distribution comparison in case of two different situations, neara hill and near a tall building. The variance in case of the hill is a lot less compared to the variance in case ofthe building thus differentiating between them. Figure 11 is a special case where the captures are taken witha corner reflector placed 2 meter from the receiver antenna and the differences in the plots are clearly visible.Table 2 provides the various moments extracted directly from the captured data. The shape parameters areextracted by fitting the empirical data to the expected distribution, thus there is a difference between thelog-normal and gamma distribution shape parameters. These correlates with the information visible in theplots and gives us a better understanding of the changes in parameters due to the change in environmentalconditions.From the plots and table provided here it is observable that there are differences in the distributions of thechannel information due to the changes in the physical environment. This enables us with the opportunity toutilise these differences and model a system that can differentiate between them. Although the differencesare quiet clear from the empirical PDF’s, it is always beneficial to further analyse the received data and gainbetter understanding.THIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272 .00 0.01 0.02 0.03 0.04 0.05 0.06 Estimated Channel Values P r o b a b ili t y D i s t r i b u t i o n PDF of north_stop_jammie_loc1 (a) Empirical PDF of captures near a bus, the receiver is locatedabout meters from the bus. Estimated Channel Values P r o b a b ili t y D i s t r i b u t i o n PDF of car (b) Empirical PDF of captures in case of single car located about meters from the receiver. Figure 8: Empirical PDF of the captured data in the presence of a bus and a car showing the differenecs inthe distributions due to the presence of particular targets.
The second set of analysis performed on the captured dataset is the Chi-square test [15, 30]. This is a testwhich takes two inputs, the observed data and the expected data and gives out a score based of the twofollowing the equation (7). The Degree of Freedom (DF) in this case is 100 since the length of the array is101 and DF is defined as one less than the length of the array. χ = M (cid:88) i =1 ( O i − E i ) E i (7)Here O i is the observed value and E i is the expected value. In order to perform this analysis PDF ofthe received data (observed data) and the distribution fitted data (expected data) are passed through equation(7) and the output of the test is recorded in Table 3. The p-value is defined as the probability of obtaining aspecific value equal to or higher than what actually is observed, maintaining the assumption that the modelis correct. The points where the p-value is null can be rejected because this means the probability of gettinga value similar or greater than that is null. All other cases in this test can be accepted and Table 3 shows thatLog-normal distribution has the most consistent result. In the next subsection we reduce the dimensionalityof the data and gain better visualization by rotating the axes. In this subsection we investigate the clustering of the estimated channel values in different scenarios. In orderto better visualize different sets of the captured data we used PCA. With this analysis we can reduce theTHIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272 .00 0.02 0.04 0.06 0.08 0.10 0.12
Estimated Channel Values P r o b a b ili t y D i s t r i b u t i o n PDF of rts_t (a) Empirical PDF of captures at a train station in the presence of atrain at a distence of 6 meters from the receiver.
Estimated Channel Values P r o b a b ili t y D i s t r i b u t i o n PDF of parking_full (b) Empirical PDF of captures at an open air parking lot, when it isfilled with cars.
Figure 9: Empirical PDF of captures at a full parking lot and near a train in a local train station showing thedifferenecs in the distributions.
Chi-Square ValueCapture Set Names Rayleigh Gaussian Log-normal Gamma
Chi-Square p-value Chi-Square p-value Chi-Square p-value Chi-Square p-value jameson stairs .
54 0 . .
61 0 . .
26 3 . e − .
97 2 . e − jameson stairs2 .
03 0 . .
34 0 . .
14 2 . e − .
09 4 . e − rts t .
03 5 . e − .
03 0 . .
29 8 . e − .
76 4 . e − roof bcr loc2 .
95 0 . .
62 0 . .
25 4 . e − .
86 5 . e − roof vh metal cr .
26 3 . e − .
82 9 . e − .
27 9 . e − .
30 3 . e − north stop jammie loc1 .
12 0 . .
97 0 . .
89 1 . e − .
07 1 . e − car .
44 8 . e − .
99 7 . e − .
07 1 . e .
73 1 . e − rhodes mem .
04 1 . e − .
23 5 . e − .
21 4 . e − .
06 1 . e − parking full .
52 1 . e − .
42 0 . .
31 2 . e − .
71 1 . e − Table 3: Chi-Square test values for the received data matched with the fitted data for each distribution.dimensionality of the data and view the datasets from an angle that provides maximum information. Here wehave observed different cluster formation of the datasets due to the change in locations and environmentalconditions. There are many different ways to reduce the dimensionality of the data and calculate the principlecomponents, for this implementation we have used Singular Value Decomposition (SVD) [3,4,27]. The majorgoals of PCA are to reduce the dimensionality of the data, extract the important information and analyse thestructure of the data. The calculations for PCA are shown in [5].In order to derive the principle components from the data we first need to generate the SVD equivalent ofthe dataset A . The input matrix A has J sets of data explained by K variables, represented by J × K . A hasa rank L with L ≤ min { J, K } , then the SVD of A will be given by equation (8). A = U ∆ V T (8) F = U ∆ (9)THIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272 .00 0.02 0.04 0.06 0.08 0.10 Estimated Channel Values -1 P r o b a b ili t y D i s t r i b u t i o n PDF of J_Stairs (a) Empirical PDF of captures taken in front of a meter tallbuilding. The receiver is located at a distance of 8 meters from thebuilding. Estimated Channel Values P r o b a b ili t y D i s t r i b u t i o n PDF of rhodes_mem (b) Empirical PDF of captures taken near a hill.
Figure 10: Empirical PDF of captures infront of a meter tall building and near a hill. Showing thedifference in the distributions between natural objects and man made objects.Here U is × L matrix of left singular vectors, V is K × L matrix of right singular vectors and ∆ is adiagonal matrix of singular values The components for PCA are obtained from the data A using equation (8).The principal component matrix F of dimension J × L is given by equation (9). To get the coefficients of thelinear combinations which are used to compute the factor scores the matrix V is used. The matrix can also beinterpreted as the projection matrix because A times V gives the projection values of the observations on theprinciple components as shown in equation (10). F = U ∆ = U ∆ V T V = AV (10)Geometrically the components can also be represented by rotating the original axes and the matrix A canbe interpreted as a product of the factor scores given by equation (11). [5] Here I is the identity matrix. A = FV T F T F = ∆ & V T V = I (11)All the plots after PCA are done by stacking the normalized datasets together and comparing the st and the nd principle components. The dataset dimensions have been reduced from 40 components to 2 major componentsand the relationship between the components are plotted to show the best view. The clusters are clearly visiblefor each data. It can be noted here that taking more components will get us more information and hence willmake scene classification even easier.Figure 13a shows the the clusters of the data captured at sea beach, near the shore, beach rock, at theTHIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272 .00 0.01 0.02 0.03 0.04 0.05 0.06 Estimated Channel Values P r o b a b ili t y D i s t r i b u t i o n PDF of roof_vh_metal_cr (a) Empirical PDF with a dihedral corner reflector. Placed in aV-H configuration, One plate was vertical and the other was kepthorizontal to the ground.
Estimated Channel Values P r o b a b ili t y D i s t r i b u t i o n PDF of roof_bcr_loc2 (b) Empirical PDF without corner reflector. Same location andclimatic conditions as Figure 11a.
Figure 11: Empirical PDF of captures with and without the presence of a dihedral corner reflector.shore but at a height of around 10 meters from the sea level. J stairs, stairs at the entrance of a building about10 meters tall. Train station, gives the captures taken at a local train station in two different situations first,when there is a train about 6 meters from the receiver and secondly in the absence of the train. Highland,refers to a location where the captures are taken on the road next to a small hilly terrain. Rhodes memorial,is situated mid way up on to a mountain in Cape Town.Figure 13b shows a particular set of captures taken with a corner reflector. The different clusters in thisfigure shows the different set of captures taken at various different configurations. CR is an acronym forcorner reflector the different allocations are in this plot. Figure 13c shows the clusters for specific set ofcaptures taken at a bus stop on the campus of University of Cape Town (UCT) at two different locations. Theuniversity shuttle service is known as ”jammie”, thus the name given. The different scenarios for this captureare, in the presence of a jammie about 4 meters away from the receiver antenna in location 1, 2 and in thesame locations when there is no jammie.Figure 13d shows the captures taken at a flat parking space when there are no cars and in the presence ofcars. This set also shows the plots for a single car dataset, where the captures are taken with and without acar in the same location on two different days (one at the morning 8 am and the other in the evening 6 pm)with different climatic conditions. Figure 12 shows the difference in the clusters at the same location whenthere is heavy rain, medium rain and on a hot day.THIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272 .0 0.8 0.6 0.4 0.2 0.0 0.2 0.4 0.6
Principle Component 1 P r i n c i p l e C o m p o n e n t PCA Index 1 & 2 comparison (Humidity, Rain and Hot Day)
Humid No RainHeavy RainHot Day
Figure 12: PCA of different climatic conditions at a single location. The conditions include high humiditybut no rain, heavy rain and Hot day with very low humidity.
Principle Component 1 P r i n c i p l e C o m p o n e n t PCA Index 1 & 2 Comparison (Beach, Highland, Building and Train)
Sea BeachBeach RockJameson StairsTrain Station No TrainTrain Station With TrainHighlandRhodes Memorial (a) The clusters shown here include the captures taken near a beach,highland, building and train.
Principle Component 1 P r i n c i p l e C o m p o n e n t PCA Index 1 & 2 Comparison (Corner Reflector (CR))
No CRNo CR loc 2Vertical CR1Vertical CR2Vertical-Horizontal CR (b) The clusters here are taken with different configuration ofCorner Reflector placed at a distance of 2 meters from the receiver.
Principle Component 1 P r i n c i p l e C o m p o n e n t PCA Index 1 & 2 Comparison (Bus)
Jammie loc2No Jammie loc2Jammie loc1No Jammie loc1 (c) Different clusters formed by captures taken near a bus indifferent configurations and locations.
Principle Component 1 P r i n c i p l e C o m p o n e n t PCA Index 1 & 2 comparison (Parking and Car)
Parking EmptyCarNo CarCar Day 2No Car Day 2Parking Full (d) Clusters of captures taken at an open air parking lot when it wasfull/empty and with a car in proximity.
Figure 13: PCA of different datasets showing the different clusters.THIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272
Conclusion
In this work we present a novel way of sensing the environment using communication signals, which wecall CommSense. The hypothesis under test was that we will be able to sense the environment by analysingthe channel equalization values of any communication system. To prove our hypothesis we chose GSMcommunication. We have implemented a very basic correlation based techniques to estimate the channelvalues from the received signal using the known training sequence. In the limited drive to prove the hypothesiswe collected data from different scenarios. On analysis it was shown that the probability distribution ofthe collected data differs as the environment changes. Hence a classic Hypothesis-test type algorithm candistinguish the scenes. Secondly, the data for different scenarios are highly clustered. This shows that a naivenearest neighbour algorithm can distinguish the environmental conditions easily. With these we prove ourhypothesis. In the future work we intend to gather more data and characterize the environment for robustclassification. We also intend to work on increasing the number of receiver nodes as another way towardsgetting robust environment sensing using CommSense.
Appendix I
In this appendix we shall briefly describe four of the most used PDF models.
The most general form of distribution used in statistical analysis is the gaussian distribution. The probabilitydensity function follows the equation (12). f ( x ) = 1 σ √ π e − ( x − µ )22 σ (12)where, µ = location parameter (mean) σ = scale parameter (standard deviation)Standard Normal Distribution where µ = 0 & σ = 1 is given as equation (13) f ( x ) = 1 √ π e − x (13)The maximum likelihood estimation for the gaussian distribution is derived by minimizing the log-likelihoodTHIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272unction of the equation (12) given by equation (14),(15) ˆ µ = 1 n N (cid:88) i =1 x i (14) ˆ σ = 1 n N (cid:88) i =1 ( x i − ˆ µ ) (15) One of the most common distributions used in wireless communication systems to model the channel isRayleigh distribution. The probability density function is given by the equation (16). f ( x ) = 1 σ xe − x σ (16)where, σ = scale parameter (mode)The maximum likelihood estimation for the parameters of rayleigh distribution is given by equation (17) ˆ σ = 12 n N (cid:88) i =1 x i x > (17) Gamma distribution is also a very versatile distribution which can be manipulated in many ways to give otherdistributions such as k-distribution which is in the scope of the future work. The general formula for theprobability density function of gamma distribution is given by equation (18) f ( x ) = ( x − µβ ) γ − exp ( − x − µβ ) β Γ( γ ) x ≥ µ ; γ, β > (18)Where, γ = shape parameter µ = location parameter β = scale parameter Γ is the gamma function given by equation (19) Γ( a ) = (cid:90) ∞ t a − e − t dt (19)THIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272he case where µ = 0 & β = 1 is called the standard gamma distribution and is given by equation (20) f ( x ) = x γ − e − x Γ( γ ) x ≥ γ > (20)The maximum likelihood estimates for the two parameter gamma distribution are calculated by solvingthe following equations (21) and (22) simultaneously ˆ β − ¯ x ˆ γ = 0 (21) log ˆ γ − ψ (ˆ γ ) − log (cid:32) ¯ x ( (cid:81) ni =1 x i ) /n (cid:33) = 0 (22)with ψ denoting the digamma function given by equation (23), which is the mathematical derivative of thegamma function. These functions are solved by using python’s scikit learn package. ψ ( z ) ≡ ddz ln Γ( z ) = Γ (cid:48) ( z )Γ( z ) (23) Distribution Shape Mean Variance Skew KurtosisGaussian
N/A . . . . Rayleigh
N/A .
25 0 .
43 0 .
63 0 . Log-normal . .
02 0 .
04 0 .
61 0 . . .
13 0 .
36 1 .
75 5 . . .
64 4 .
67 6 .
18 110 . . .
39 2926 .
35 414 .
36 9220556 . Gamma . . . .
83 12 . . . . .
41 3 . . . . .
89 1 . . . . .
67 0 . Table 4: Moments of different distributions shown in Figure 14. If x is a random variable distributed lognormally then y = ln( x ) is normally distributed where ln in naturallog. It is very helpful in distinguishing the dataset in this work because it gives the best fit so far. Theprobability density function is given by equation (24) f ( x ) = e − ((ln(( x − θ ) /m )) / (2 σ )) ( x − θ ) σ √ π x > θ ; m, σ > (24)THIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272 -40 -37 -34 -31 -28 -25 -22 -19 -16 -13 -10 -7 -4 -1 P r o b a b ili t y Gaussian Distribution µ =0 ,σ =0 . µ =0 ,σ =1 . µ =0 ,σ =5 . µ =1 ,σ =1 . µ = − ,σ =1 . -86 -79 -72 -65 -58 -51 -44 -37 -30 -23 -16 -9 -2 P r o b a b ili t y Rayleigh Distribution σ =0 . σ =1 . σ =2 . σ =3 . σ =4 . X -17 -15 -13 -11 -9 -7 -5 -3 -1 P r o b a b ili t y Gamma Distribution k =0 . , Θ=1 . k =1 . , Θ=2 . k =2 . , Θ=2 . k =3 . , Θ=2 . k =5 . , Θ=1 . k =7 . , Θ=1 . k =9 . , Θ=0 . X -6 -5 -4 -3 -2 -1 P r o b a b ili t y Lognormal Distribution µ =0 . ,σ =0 . µ =0 . ,σ =0 . µ =0 . ,σ =1 . µ =1 . ,σ =2 . µ =2 . ,σ =5 . Figure 14: All distributions explained above (y-axis is in log scale).where, σ = shape parameter (standard deviation) θ = location parameter m = scale parameter (median)The standard lognormal distribution is given when θ = 0 & m = 1 and is denoted as equation (25) f ( x ) = e − ((ln x ) / σ ) xσ √ π x > σ > (25)The lognormal distribution is commonly characterised with its mean µ given as µ = log m , using thisparameter the density function changes to equation (26) f ( x ) = e − (ln( x − θ ) − µ ) / (2 σ ) ( x − θ ) σ √ π x > σ > (26)THIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272aximum likelihood estimation of this distribution is given by equation (27),(28),(29) ˆ µ = 1 N N (cid:88) i =1 ln X i (27) ˆ σ = 1 N N (cid:88) i =1 (ln ( X i ) − ˆ µ ) (28) ˆ m = exp ˆ µ (29)Table 4 and Figure 14 shows the moments of the above mentioned distribution. Appendix II
The capture named “J Stairs” or “jameson stairs” or “jameson stairs2” is taken at a location in UCT as shownin Figure 15.
J stairs capture taken at this location
Figure 15: Location for capture named “J Stairs” or “jameson stairs” or “jameson stairs2”THIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.150272 eferences [1] Digital cellular telecommunications system (phase 2+); multiplexing and multiple access on the radiopath. ETSI TS 145 002 V12.3.0, (2015-01). (3GPP TS 45.002 version 12.3.0 Release 12).[2] Digital cellular telecommunications system (phase 2+); radio transmission and reception. ETSI TS 145005 V12.4.0, (2015-01). (3GPP TS 45.005 version 12.4.0 Release 12).[3] Herv´e Abdi. The eigen-decomposition: Eigenvalues and eigenvectors.
Encyclopedia of measurementand statistics , pages 304–308, 2007.[4] Herv´e Abdi. Singular value decomposition (svd) and generalized singular value decomposition.
Encyclopedia of measurement and statistics. Thousand Oaks (CA): Sage , pages 907–912, 2007.[5] Herv´e Abdi and Lynne J Williams. Principal component analysis.
Wiley Interdisciplinary Reviews:Computational Statistics , 2(4):433–459, 2010.[6] Christian R Berger, Bruno Demissie, J¨org Heckenbach, Peter Willett, and Shengli Zhou. Signalprocessing for passive radar using ofdm waveforms.
Selected Topics in Signal Processing, IEEE Journalof , 4(1):226–238, 2010.[7] C Coleman and Heath Yardley. Passive bistatic radar based on target illuminations by digital audiobroadcasting.
Radar, Sonar & Navigation, IET , 2(5):366–375, 2008.[8] GNU Radio Website. Gnu radio website. , accessed October 2015.[9] Hugh Griffiths. Passive bistatic radar.
RTO Educational Notes-Lecture Series RTO-EN-SET-133Waveform Diversity for Advanced Radar Systems, Brno , 2009.[10] Friedhelm Hillebrand.
GSM and UMTS: the creation of global mobile communication . John Wiley &Sons, Inc., 2002.[11] J Homer, K Kubik, B Mojarrabi, ID Longstaff, E Donskoi, and M Cherniakov. Passive bistatic radarsensing with leos based transmitters. In
Geoscience and Remote Sensing Symposium, 2002. IGARSS’02.2002 IEEE International , volume 1, pages 438–440. IEEE, 2002.[12] Michael Inggs, Cunsheng Tong, Roaldje Nadjiasngar, Graham Lange, Anadi Mishra, and FrancoisMaasdorp. Planning and design phases of a commensal radar system in the fm broadcast band.
Aerospace and Electronic Systems Magazine, IEEE , 29(7):50–63, 2014.THIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.15027213] Michael Inggs, Andrew van der Byl, and Cunsheng Tong. Commensal radar: Range-doppler processingusing a recursive dft. In
Radar (Radar), 2013 International Conference on , pages 292–297. IEEE, 2013.[14] Michael R Inggs, CA Tong, Akhilesh Kumar Mishra, and FDV Maasdorp. Modelling and simulationin commensal radar system design. In
Radar Systems (Radar 2012), IET International Conference on ,pages 1–5. IET, 2012.[15] MG Kendall and A Stuart. The advanced theory of statistics vol. 2, inference and relationship. charlesgriffin and co., ltd, 1961.[16] P Krysik, K Kulpa, M Baczyk, Ł Ma´slikowski, and P Samczynski. Ground moving vehicles velocitymonitoring using a gsm based passive bistatic radar. In
Radar (Radar), 2011 IEEE CIE InternationalConference on , volume 1, pages 781–784. IEEE, 2011.[17] Amit Kumar Mishra. Monitoring changes in an environment by means of communication devices,April 20. PCT Application No: PCT/IB2016/052235.[18] Amit Kumar Mishra. Commsense: Radar system using existing communication infrastructure to sensethe environment. In
Radioelektronika (RADIOELEKTRONIKA), 2015 25th International Conference ,pages 276–279. IEEE, 2015.[19] Amit Kumar Mishra and Santu Sardar. Application specific instrumentation and its feasibility foruwb sensor based breast cancer diagnosis. In
Power, Control and Embedded Systems (ICPCES), 2010International Conference on , pages 1–4. IEEE, 2010.[20] Nuand Website. Nuand website. http://nuand.com/ , accessed October 2015.[21] Liang Pu, Jian Liu, Yuan Fang, Wei Li, and Zhisen Wang. Channel estimation in mobile wirelesscommunication. In
Communications and Mobile Computing (CMC), 2010 International Conferenceon , volume 2, pages 77–80. IEEE, 2010.[22] Markku Pukkila. Channel estimation modeling.
Nokia Research Center , 2000.[23] P Samczynski, K Kulpa, M Malanowski, P Krysik, and L Ma´slikowski. A concept of gsm-based passiveradar for vehicle traffic monitoring. In
Microwaves, Radar and Remote Sensing Symposium (MRRS),2011 , pages 271–274. IEEE, 2011.[24] Santu Sardar and Akhilesh Kumar Mishra. ASIN-based UWB radar for sludge monitoring.
Access,IEEE , 2:290–300, 2014.THIS PAPER IS PUBLISHED IN IEEE AES MAGAZINE, VOLUME 32 ISSUE 2,UNDER DOI: 10.1109/MAES.2017.15027225] Santu Sardar and Amit Kumar Mishra. LTE-CommSense: LTE communication infrastructure basedsensing for environment monitoring. In . IEEE, 2015.[26] Hongbo Sun, Danny KP Tan, and Yilong Lu. Aircraft target measurements using a gsm-based passiveradar. In
Radar Conference, 2008. RADAR’08. IEEE , pages 1–6. IEEE, 2008.[27] Yoshio Takane. Relationships among various kinds of eigenvalue and singular value decompositions.In
New developments in psychometrics , pages 45–56. Springer, 2003.[28] Danny KP Tan, Hongbo Sun, Yilong Lu, M Lesturgie, and HL Chan. Passive radar using global systemfor mobile communication signal: theory, implementation and measurements. In
Radar, Sonar andNavigation, IEE Proceedings- , volume 152, pages 116–123. IET, 2005.[29] Danny KP Tan, Hongbo Sun, Yilong Lu, and Weixian Liu. Feasibility analysis of gsm signal for passiveradar. In
Radar Conference, 2003. Proceedings of the 2003 IEEE , pages 425–430. IEEE, 2003.[30] Ali Vardasbi, Mahmoud Salmasizadeh, and Javad Mohajeri. Multiple-chi-square tests and theirapplication on distinguishing attacks. In