Crowd monitoring systems enhance the security of large public events. Event managers have expressed their interest in leveraging modern technologies to i) monitor events in real time [1, Sec. 7] , ii) predict crowd counts in the future [1, Sec. 5.1.1], and iii) perform post-analyses. In particular, computing real-time crowd densities in strategic areas allows security managers to decide whether an event has reached its maximum capacity [1, 2] . Crowd count time series can be fed into forecasting algorithms to predict overcrowding-which allows security personnel to execute countermeasures anticipatedly. Our previous works [3, 4] deal with the issue of forecasting crowd counts.
Such systems also have other use cases, albeit similar to those of the previous paragraph. Within the framework of our partnership with Brussels Major Events-a structure responsible for public events in Brussels-we installed a crowd monitoring system on one of the main commercial streets of Brussels, namely Rue Neuve (Nieuwestraat in Dutch), to estimate attendance during winter sales. It has been reinstalled in the same street to track attendance as Covid-19 lockdown measures get incrementally relaxed. Finally, we also installed our monitoring system in the largest library of our university: the Humanities library. Generally, our monitoring system can be used to monitor indoor and outdoor public places, if they are sufficiently large and
A recent work is [5] (and the related one [6] ). The authors deployed tens of nodes across rooms to be monitored and make them communicate with one another. The received powers for all communication links are then a proxy for the number of attendees, because human bodies attenuate WiFi signals. This solution is fully non-cooperative, is compatible with low numbers of attendees (¡100 people), is not affected by MAC address randomization and can be calibrated easily when the monitored room is empty. However, nodes must be at a low height (¡ 2 meters) for human bodies to attenuate signals. Moreover, tens of nodes are necessary to monitor a single room (we estimated they installed one node per 15 to 40 square meters based on [5, Fig. 2, 13, and 24] ). Besides, their counting errors are higher than ours (see, e.g., [5, Fig. 11] and other figures in their paper). Similarly, some researchers also rely on channel state information (CSI) to identify human behavior [7] . Other seminal and recent works relying on CSI are [8] [9] [10] [11] .
The work [12] proposes a crowd monitoring solution for user localization in large buildings. They rely on clients connected to access points to which they have access. Therefore, their approach is partially cooperative. As a result, they depend on users using their access points but do not have to deal with MAC randomization issues. A similar work is [13] .
Another work, [14] , presents a WiFi sensing system for monitoring attendance and travel times of a bus line in Charlottesville, Virginia. It seems the estimator in [14] discards MAC addresses appearing only once [14, Sec. 3.2.1] , thereby filtering out randomized MAC addresses (which were likely less frequent in 2017 [14, Fig. 4 ] than in 2019, the year we collected the measurements presented in this paper). Their solution has been tested for low numbers of people whereas ours is suited to a larger number of attendees. [15] and [16] are similar works that describe a system tracking buses using WiFi signals emitted by its riders to predict their arrival times.
A work similar to ours is [17] , whose described system collects data essentially identical to ours (entries that consist of a timestamp, a MAC address and a received signal strength indicator). Their focus is on density monitoring and trajectory tracking. They do not refer to MAC address randomization, probably because their measurements were obtained a few years ago (between 2014 and 2016 according to [17, Sec . IV-C]), a time for which MAC address randomization was not a significant issue. Therefore, it is not clear that the accuracy of their monitoring system would be as high with today's smartphone anonymization. Our system is not affected by MAC address randomization and we provide a privacy analysis of our system as well as a theoretical analysis of its accuracy.
An older and seminal work is [18] in which the authors emulate APs for common service set identifiers (SSIDs) and SSIDs present in the information elements (IEs) of detected PRs. They also send request to send (RTS) packet injection. [18] thus describes an active scanning system. In [19] , the authors propose an indoor localization system using jointly micro-electro-mechanical systems and WiFi fingerprinting.
The authors of [20] and [21] presented a crowd monitoring based on WiFi probe requests. Key differences with our work are that i) our estimator uses all detected probe requests whereas that of [21] filters out all locally administered MAC addresses [21, Sec. 4] , ii) our system anonymizes MAC addresses in such a way that tracking is not possible (at least for more than one minute), in order to be out of GDPR, iii) our system uses SHA-256 hashing in conjunction with cryptographic peppers while the other work uses SHA-256 without peppers [21, Sec. 3.3] , thereby making their anonymization procedure more vulnerable to brute force attacks [22] [23] [24] , and iv) their system improves localization accuracy using Friis' formula [21, Sec. 4.5.1] , relying on the assumption that the distance is statistically correlated with the received power in an omnidirectional way; it is a good approximation only if line-of-sight (LOS) paths are dominant.
An excellent review of recent works in crowd monitoring making use of WiFi is [21, Sec. 2 and Table 1 ]. Other more general reviews are [25] , [12, Sec. 1.1] , [14, Sec. 2 and Table 1 ] and [14, Sec. 2] . We also wrote a more detailed review of slightly older works in [3, Sec. I.B], which we do not rewrite here for the sake of brevity.
We gave a minimal overview of this system in our previous works on forecasting [3, 4] , whose main purpose was to demonstrate the interest of crowd monitoring systems for forecasting. We also compared our counts against those of a telco operator in ideal conditions to show our extrapolation factor is realistic. This manuscript details more deeply the system architecture and compares our counts against those from a camera-based system for an indoor environment. It also presents mathematical results on the accuracy of the estimator and the effects of our anonymization procedure. Finally, it presents a detailed complexity analysis of our counting algorithm.
The notations P and E denote the probability of an event and the mathematical expectation, respectively. Lower-case bold letters denote (column) vectors and upper-case bold letters denote matrices. The inner product of x and y is x, y .
The estimated crowd counts that we use for forecasting purposes are derived from PRs [26, Chapter 4] . WiFi devices periodically transmit PRs to request nearby access points (APs) to send back probe responses. This is an active scanning mechanism to discover APs. WiFi devices transmit PRs even when not linked to a WiFi network. Thus, measuring the rate of PRB transmission in an area gives an idea of the number of WiFi-enabled devices in the covered area, a number which can be extrapolated to a crowd count. See Figure 1 for an illustration of the process. Several almost identical PRs are sent in a row, within a time frame lasting less than 10ms [27, Sec. 2.1]; in this paper, we call a set of such PRs a probe request burst (PRB). 
PRs contain a source address (SA) field of six bytes [26, , which is usually a randomized MAC address. Recent and updated operating systems indeed embed this randomization process to make smartphone tracking difficult [27] [28] [29] .
Some older works from 2016-2017 show that anonymized PRs may be "deanonymized" (see, e.g., [29, 30] ). In the future, however, deanonymization methods may not work if operating systems strengthen anonymization. For example, [29, Section 4] partially relies on sequence numbers [26, Figure 4 -52], which are numbers associated with each PR that are incremented in between consecutive PRs. So far, it appears that such sequence numbers are not randomly reset from one PRB to the next one-a fact that the authors [26] leverage to track smartphones. Should sequence numbers be randomly reset in the future, the strategy may not work anymore. More generally, MAC address randomization is likely to get tougher in the future [5, Sec. 1]; as pointed out in [21, Sec. 4] , "the IEEE 802.11 working group has created a Topic Interest Group (TIG) on Randomized and Changing MAC addresses (RCM)".
As discussed later on, MAC address randomization does not affect our counts, which makes our system future-proof, in opposition to other WiFi monitoring systems either deanonymizing PRBs or identifying non-randomized PRBs (see [21] ).
We assume the monitored area hosts n ppl individuals. Each individual has a probability p i ≤ 1 (1 ≤ i ≤ n ppl ) to transmit a PRB within an elementary time frame, which is zero if the individual's smartphone has no WiFi enabled. This probability exists because, for a fixed smartphone configuration, the initial time at which a smartphone starts transmitting PRBs is random and follows a uniform distribution. For example, for time frames with a duration of 10 seconds and for a smartphone transmitting PRBs every 30 seconds, p i = 1/3. Equivalently, p i may be thought of as a rate of PRB transmission divided by the duration of an elementary frame.
Each probability p i is fixed by the associated smartphone configuration (whether WiFi is enabled and how frequently it sends PRBs). Individuals' smartphones (and their configurations) can be considered to be independently and identically distributed (iid) random variables, which makes the {p i } i iid random variables. We have K < ∞ possible values {α k } 1≤k≤K for p i ; the probability r k := P[p i = α k ] obeys K k=1 r k = 1, with α k = 0 corresponding to having no WiFi-enabled device.
The number of distinct PRBs within a time frame is X :=
We propose an unbiased estimator of the number of individuals that isĈ
where E[Ĉ] = n ppl with extrapolation factor β :
The variable X is a scaled sum of n ppl statistically independent and identically distributed (iid) Bernoulli random variables X i of parameter p := E[p i ]. As a result,Ĉ/β follows a binomial distribution B(n ppl , p).
We now present a concentration inequality for our estimatorĈ around its mean. Loosely speaking, this inequality is theoretical evidence our estimator is reliable. We use results derived in [31] and compare them against a canonical concentration inequality for bounded random variables. A key quantity depending on p is K(p), defined below.
with q := 1 − p.
Proposition 1 helps understanding the shape of K(p).
With K defined as in (2), we have the following properties: Let us now state the concentration inequality forĈ.
Proposition 2. With K defined by (2) andĈ by (1), we have, for any ϕ > 0,
Proof. We have already established that X is of a sum of n ppl iid. Bernoulli random variables of parameter p. Thus, [31, Corollary 6.1 (ii)] directly implies
.
This concentration inequality upper bounds the probability thatĈ diverges from its mean n ppl as a function of a proportion ϕ of the mean. In particular, it shows that the probability of a divergence of ϕn ppl decreases exponentially with the number of people in the area n ppl -which means that the relative accuracy of our estimator increases with n ppl and becomes infinite as n ppl → ∞.
As lim p→1 − p 2 K(p) −1 = +∞, if every individual is guaranteed to send one PR during any elementary time frame, our relative estimator accuracy is infinite. Conversely, using L'Hôpital's rule, we also have
which suggests that with no individual ever sending a PRB, the estimator is worthless.
We now compare our bound against Hoeffding's inequality (see [32] and [33, Theorem 2.8] ). Without proof details, we easily obtain Hoeffding's inequality:
which is also obtained by using (3) and K(p) ≤ 1/4 (see [31, Remark 5 .1]), which shows (3) outperforms (4).
Our system comprises i) a set of sensors, ii) a processing subsystem on a central server collecting all PRBs and processing them in real time, and iii) a dumping subsystem (that is part of the central server) that further anonymizes and then dumps PRBs. All communications between the sensors and the central server use layers of authentication; they are secured using HTTPS, thereby encrypting packets and also preventing man-in-the-middle attacks.
We developed sensors i) detecting PRBs, ii) anonymizing them and iii) sending them to a central server.
• A Raspberry Pi 3B (with Raspbian Stretch).
• An Alfa AWUS036NHA WiFi dongle (chipset Atheros AR9271L) supporting monitor mode-a state that makes the dongle capture all over-the-air WiFi messages, without being restricted to those of a particular WiFi network. We use the dipole antennas shipped with Alfa AWUS036NHA dongles. Sensor antennas point perpendicularly to the ground.
• A 4G dongle granting access to the Internet Figure 2 shows a photograph of our sensor.
We wrote our sniffing program in C++, which is multi-threaded and uses packet capture library libpcap. For each detected PRB, our sensors send [3, Sec. II-B] "i) an anonymized MAC address of the PRB, ii) the timestamp of detection iii) a received signal strength indicator (RSSI) value, which is a number quantifying the received power ". We stress tested our sensors to ensure neither the WiFi dongle nor the Raspberry Pi fail to handle large PRBs transmission rates. 
All sensors periodically retrieve from the central server an up-to-date array of (cryptographic) server peppers. Each pepper of the array is associated with a one-minute time frame, during which it will be used. The central server regenerates the server peppers in real time, and it deletes old peppers so that they cannot be retrieved in the future. The server uses an entropy pool (/dev/urandom on Linux distributions) to generate cryptographically secure peppers. A sensor pepper is also hardcoded in the C++ codebase of all sensors; it is common to all sensors (at least all sensors located in the same area and thus likely to detect identical PRBs simultaneously). It is a final line of defense in case the server peppers get compromised.
As depicted in Figure 3 , for every received probe request, the sensor prepends a global pepper to the full MAC address before computing the SHA-256 hash of the concatenated byte sequence, whose 256 bits are truncated to 64 bits. The pepper is the concatenation of the sensor pepper and the server pepper, both of 128 bits.
As shown in [34] , the system we designed for anonymization purposes meets four essential requirements [34, Sec. II]. First, time synchronization is accurate enough to make sensors use identical peppers at identical time instants (at least when operating on networks with low latency, such as LTE networks [35] ). Second, from the SA identifiers, it is realistically impossible to recover the original MAC addresses [34, Sec. II.C]. Third, tracking individuals for more than one minute is not possible [34, Sec. II.D]. Fourth, the collision rate of the truncated SHA-256 hash is less than 10 −9 for 10 million MAC addresses (which corresponds to an unrealistically high number of individuals). Satisfying the first and fourth requirements ensures anonymization does not tamper with the counting method. The second and fourth requirements consist in privacy-enhancing features. 
For any elementary time frame T k of one minute, we extract from our database all the anonymized MAC addresses whose timestamps belong to [t k − T + 1, t k ] (with T = 60 s), which creates an array arr mac of 3-tuples; the ith 3-tuple is (S (i) , aMAC (i) , RSSI (i) )-where S (i) is a sensor ID, aMAC (i) denotes the ith anonymized MAC address and RSSI (i) is the ith RSSI. Each anonymized MAC address may appear several times because sensor ranges usually overlap.
The first step consists in computing counts for every sensor. We initialize a hash table whose keys and values are anonymized MAC addresses and 2-tuples (sensor ID, RSSI), respectively. We go through all the rows of arr mac.
Before taking into account a row from arr mac, we apply a sensor-specific RSSI thresholdwhich, if violated, discards the row. This consists in a lower bound on the RSSI that artificially reduces the range of a sensor. We used such thresholds in our experiments, see Sec. IV. We then check if the anonymized MAC address of the row exists in the hash table. If it exists, and if the stored RSSI is less than that of the current row of arr mac, we replace the 2-tuple value in the hash table with a new 2-tuple corresponding to the sensor ID and the RSSI of the row.
After having looped through arr mac, the hash table contains, for each anonymized MAC address, the sensor ID with the highest observed value of the RSSI. Thus, a count can be generated for each sensor and counts of different sensors take into account any MAC address only once. The count of any area is then computed by summing the counts of the sensors that area indexes. Algorithm 1 describes the procedure.
Require: List of PRBs arr mac for the time frame of interest only, 3-tuples (sensorid, amac, rssi);
List of sensors sensors, 2-tuples (sensorid, rssilowerbound) 1: Initialize a hash table ht whose keys are 64-bit anonymized MAC addresses (or equivalently, random tokens) and whose values are 2-tuples (sensorid, rssi 
We now turn to the complexity analysis (in time and in memory). We assume that there are n S sensors, each one of which captures no more than n meas PRBs for a time frame. Storing all the PRBs for a given time frame has a memory footprint of n S n meas 16 10 −6 MB.
The memory of the hash table used for processing PRBs is also reasonable. Let n b denote the number of buckets of the hash table. In practice, n b can be chosen to get a load factor lower than or equal to α. Thus, let us fix n b = n S n meas α −1 . Setting the number of buckets beforehand requires us to know approximately the maximum number of PRBs per time frame attained in practice (and the proportion of duplicated PRBs). Rehashing as the hash table grows is an alternative. In practice, because the SA identifiers are approximately uniformly distributed, no Figure 4 : Example of a C structure representing a probe request burst. In this case, any standard compiler appends 1 trailing pad byte for data alignment purposes; thus, the structure size is 16 bytes. The size of rssi is that of the antenna signal field of standard RadioTap headers. hash function is needed and the rehashing has a lower cost than in a textbook hash table.
With a C structure similar to that of Figure 4 , each 2-tuple (sensor ID, RSSI) of the hash table consists of 8 bytes (including two trailing pad bytes). Assuming that collision resolution relies on separate chaining with linked lists [36, Chap. 11] , the baseline memory footprint of the hash table is equal to n b 8 10 −6 MB on a 64-bit architecture. Every node of the linked list has a memory footprint of 16 bytes (8 bytes for the pointer and 8 bytes for the 2-tuple value). Thus, when loading n S n meas α −1 entries in the hash table, its memory footprint is of n S n meas (8 α −1 + 16) 10 −6 MB (the first and second terms correspond to the bucket pointers and the nodes of the linked lists, respectively).
As a conclusion, processing a large event is computationally tractable from a memory complexity point of view. For large events lasting several days, loading all the measurements in memory at once may be impossible but is also pointless; the proposed method processes time frames T k sequentially and independently from one another.
With a properly designed hash table, insert and search operations have an average time complexity of O(1). Looping through all entries in arr mac has a time complexity of O(n S n meas ). The reason is that the number of loops is n S n meas , each one of which includes a search operation, a comparison, and possibly an insert or a replace operation (all these operations are O (1)). Counting all entries in the final hash table with specific sensor IDs has a time complexity of O(n S n meas ) because the prescribed load factor makes the number of buckets directly proportional to n S n meas . Releasing the linked lists of all buckets also has a time complexity of O(n S n meas ); the program performs this operation in between time frames to reset the hash table. Globally, the average time complexity is O(n S n meas ). It is easy to show that the worstcase time complexity is O((n S n meas ) 2 )-for the extreme case for which all the SA identifiers are mapped onto the same bucket, thereby creating a unique linked list of size n S n meas .
Our system periodically dumps PRBs from the SQL database into binary files. Each dump file corresponds to a particular sensor and a particular day. This allows us to keep in check the size of the SQL table storing PRBs and its indexes. It also enables us to store these files in a cheap storage location (e.g., in "cold storage" facilities) or on on-premises drives for backup.
The SQL database stores anonymized MAC addresses; theoretically, a deterministic link still exists between the original MAC address and its corresponding SA identifier. Removing the link is beneficial because someone could identify a vulnerability of SHA256 in the future. Therefore, our dumping program randomizes SA identifiers per time frame T k using, e.g., the Mersenne twister. The links "SA identifier → final SA identifier" is reset after each time frame. A cryptographically secure pseudorandom number generator (CSPRNG) is not needed as our only requirements are i) to remove any deterministic link between the original MAC address and the identifier ii) and having uniformly distributed identifiers. This approach also makes it impossible for hackers to revert their way back to the original SAs on the basis of the dump files, even if they intercept the peppers.
A previous evaluation focusing on the extrapolation factor for public events is [3, Sec. II-E and Fig. 2] ; this former analysis compares our WiFi counts with those from a telco operator. This paper and section provide an experimental evaluation of the accuracy of our WiFi system in an indoor environment this time. We rely on third-party counts from Affluences and their 3D Video sensor system [37] , which has been installed at the entries and exits of the Humanities library at Université libre de Bruxelles (ULB).
As in [3] , we use two accuracy measures: the root mean square error (RMSE) and the mean absolute percentage error (MAPE). For a time series {x t } 0≤t≤N −1 of N true counts and a time series {x t } 0≤t≤N −1 of N approximated counts,
and
Our setup at the Humanities library consists in six sensors installed on three (consecutive) floors of an eight-story building.
In ideal circumstances, sensors cover the whole area to be monitored. In practice, budget or infrastructure constraints may prevent an installation with a full coverage and the total counts of people is extrapolated on the basis of counts for a sub-area. Thereby, withĈ (part) denoting the (partial) counts for the covered sub-area, we havê
whereĈ and β are defined in Section II.C and κ is an extrapolation factor converting counts for the sub-area into counts for the whole area. (If the whole area is covered, κ = 1.) The global extrapolation factor is thenβ := κβ. A more complete model that includes noise signals for both extrapolations isĈ
where C denotes the true count whereas e (κ) and (β) denote errors linked to the two extrapolation procedures. As the cameras provide counts for the whole library and the WiFi system covers three stories of out eight, we have κ > 1 and e (κ) = 0.
We shall fit the global extrapolation factorβ using a least squares approach with N measurements for each subsystem. Let c Affl. ∈ R N and c WiFi ∈ R N denote counts from the Affluences cameras and WiFi subsystems, respectively. Affluences counts are available every 30 minutes and we compute WiFi counts every 5 minutes. We downsample the WiFi count series by 6 to obtain comparable and compatible time series for both subsystems. We particularize the linear model y = Ax using the substitutions y = c Affl. , A = c WiFi and x = [β] ∈ R 1×1 . The pseudoinverse of A with linearly independent columns is
which provides a least squares estimate forβ that is
IV.D Preprocessing pipeline
Our WiFi system works all the time; however, its accuracy should only be evaluated during opening times. To do this, we preprocess both the Affluences and WiFi time series in the following way: 3. Restrict the time ranges to those during which the library is guaranteed to be opened (from 9:00 AM to 6:00 PM).
First of all, we make the pessimistic assumption thatβ is constant over time. This is not necessarily true because we only cover three floors of the library and students pursue different endeavors over time; for example, many projects are over by May, which means that the students spread differently in the floors of the library as they use less frequently the rooms to discuss with fellow classmates for projects. This pessimistic approach gives us a lower bound on the accuracy of the WiFi system because e (κ) = 0 and e (κ) is an error term linked to the partial coverage that does not usually appear for ideal installations. In other words, any system with limited coverage would be subject to noise e (κ) and all systems with full coverage have e (κ) = 0. Figure 5 compares Affluences counts against WiFi ones, for the restricted time frame running from 09:00 AM to 6:00 PM. Figure 6 does the same but displays the full days, which makes the plot easier to read.
The estimated global extrapolation factor is equal to 5.031 (for time frames of T = 60 seconds), see (8) . Comparatively, in our previous studies with full coverage [3, 4] , we obtained a value of 3 for time frames of T = 30 seconds (which is equivalent to an extrapolation factor of 1.5 with T = 60 seconds). This suggests κ 5/1.5 3.33, which is realistic given our coverage (three floors out of eight). For Figure 5 , the RMSE and MAPE values are of 120.9 and 12.7 %, respectively. The mean of the counts is equal to 824. These figures are thus upper bounds on the error of the WiFi system. Globally, our accuracy estimate based on indoor measurements is pessimistic for large events or buildings because the experiment we could carry out suffers from errors linked to:
1. the relatively low number of people in the monitored area (300 people on the three stories against thousands in larger events)
2. the extrapolation of the crowd counts on three floors to eight floors 3. the use of RSSI thresholds that we have been tuned coarsely. In large events or using directional WiFi antennas, however, such thresholds would not be necessary as the whole area is large and surrounding areas do not host a significant number of attendees.
We now estimate the global extrapolation factor for each week separately, to better reflect the time-varying distribution of the students across the different floors. Mathematically, it translates into a partial extrapolation factor κ in (7) being a function of time. We point out again that the time-varying nature of the extrapolation factor stems purely from ours monitoring a sub-area and extrapolating counts to the full area. The global extrapolation factor β is constant for events or buildings that are fully covered.
Albeit a rather theoretical exercise, compensating the time-varying nature of the extrapolation factor gets our accuracy figures closer to those that would have been obtained with full coverage. The improvements resulting from this exercise also suggest that our partial coverage leads to inflated errors in comparison to full-coverage scenarios. 
Measured counts -WiFi Measured counts -Affluences Figure 6 : Full comparison of camera and WiFi counts on selected days, with a global extrapolation factor of 5.031. Table 1 reports the results, including the ones of Sec. IV.E on its last row. We observe that the global extrapolation factor estimates increase over time, which stems from the humanities library becoming more crowded as examination sessions get closer. A possible explanation is that students favor working on floors that happen to be covered by the sensors and move to the remaining floors as seating options become scarcer; thus, the global extrapolation factor increases over time.
Finally, as expected, using extrapolation factors optimized per week improves the average RMSE and MAPE in comparison to using the one obtained for the global, aggregated time series. Nevertheless, the RMSE and MAPE improvements stemming from using weekly extrapolation factors is lower than 10 %. 
In public events, our experience is that sensors are often not connected to a dedicated power supply line, sharing instead power supplies with other devices (e.g., lightning devices). These other circuits may be unplugged to save power at night or during daytime. Even if the sensors were connected to dedicated circuits, these could malfunction or be shut down for maintenance without prior notice. Therefore, we recommend making sensors unaffected by improper shutdowns, by using high-quality persistent storage (e.g., using high-end eMMC memory) and by mounting the operating system in read-only mode.
This paper describes a crowd monitoring system relying on probe requests transmitted by attendees' smartphones in the monitored area. This system is suitable for indoor and outdoor areas hosting at least a few hundred attendees. The monitoring system ensures strict privacy requirements are met and is therefore compatible with modern privacy laws. We provided both theoretical and experimental evidence that our system computes accurate estimates of the number of attendees. Despite non-ideal experimental conditions, the MAPE we computed is of less than 13 %.