Determining locations of mobile devices from wireless signals BACKGROUND OF THE INVENTION This invention relates to methods, software and apparatus for determining locations of mobile devices of a network of devices. It is known to determine the locations of mobile devices, such as RFID tags or cell phones, of a network of devices (which may also include one or more static devices), by transmitting acoustic or radio signals between the devices in order to estimate the distances between pairs of the devices. Locations may be estimated using techniques such as trilateration or pattern-matching. However, it is not straightforward to do so efficiently and reliably. The present invention seeks to provide a novel approach for determining the locations of mobile devices efficiently and reliably. SUMMARY OF THE INVENTION From a first aspect, the invention provides a method of determining a respective location of each of a set of one or more mobile devices of a network of devices, the method comprising: receiving data comprising a set of values representative of distances between respective pairs of devices of the network of devices, determined from wireless signals transmitted between the devices; processing the values representative of distances to determine a set of coefficients for a system of linear equations; and solving the system of linear equations to determine two-dimensional or three- dimensional coordinates representative of the location of each of the set of one or more mobile devices in a two-dimensional plane or a three-dimensional space. From a second aspect, the invention provides computer software for determining a respective location of each of a set of one or more mobile devices of a network of devices, wherein the computer software comprises instructions which, when executed by a processing system, causes the processing system to: receive data comprising a set of values representative of distances between respective pairs of devices of the network of devices, determined from wireless signals transmitted between the devices; process the values representative of distances to determine a set of coefficients for a system of linear equations; and solve the system of linear equations to determine two-dimensional or three- dimensional coordinates representative of the location of each of the set of one or more mobile devices in a two-dimensional plane or a three-dimensional space. From a third aspect, the invention provides a location system comprising a processing system and a memory, wherein the memory stores software comprising instructions which, when executed by the processing system, cause the processing system to: receive data comprising a set of values representative of distances between respective pairs of devices of the network, determined from wireless signals transmitted between the devices; process the values representative of distances to determine a set of coefficients for a system of linear equations; and solve the system of linear equations to determine two-dimensional or three- dimensional coordinates representative of the location of each of a set of one or more mobile devices in a two-dimensional plane or a three-dimensional space. The processing system may, in some embodiments, further comprise the network of devices, which may include one of more static devices in addition to the one or more mobile devices. Thus it will be seen that, in accordance with embodiments of the invention, two- dimensional (2D) or three-dimensional (3D) coordinates can be determined for one or more mobile devices by solving a linear system. This approach is based on the insight that the locations of the one or more devices, considered as coordinates in the plane or in a three-dimensional space, can be represented by a consistent set of linear combinations of the locations of the other devices in the network, leading to a linear relationship between the device locations. This linearity enables a linear system to be solved to determine the locations of individual or multiple devices efficiently and consistently. It can also advantageously enable the locations of mobile devices to be determined accurately even for devices that are in effective wireless range only of other devices that are also of unknown location—i.e. even when there is no direct distance measurement available from the mobile device to a device of fixed or known location, such as to a static beacon or a mobile device of known position. It will be appreciated that, in some embodiments, the coordinates may be two- dimensional coordinates (e.g. x, y) representative of a respective location of each of the one or more mobile devices in a two-dimensional plane (e.g. corresponding to a plan view of an environment). However, in some embodiments, the coordinates may be three-dimensional coordinates (e.g. x, y, z) representative of respective locations of the one or more mobile devices in three-dimensional space (e.g. including height information). In some embodiments, the wireless signals may be radio signals, such as Bluetooth Low Energy™, WiFi™ or NFC signals. In some embodiments, the wireless signals may be acoustic signals, such as ultrasound signals. A combination of radio and acoustic signals may be used to determine the data representative of distances. In some embodiments, one or more of the mobile devices may be a smartphone or an electronic tag configured to send and receive wireless signals. The network of devices may be arranged to communicate as a radio mesh network. In some embodiments, data representative of distances between pairs of devices of the network may be determined based on a signal strength measure of the wireless signals transmitted between the devices. The data may be determined based on one or more of received signal strength indicators (RSSI) of the wireless signals, angle of arrival (AoA) measurements, round trip time (RTT) measurements, time of arrival (ToA) measurements, or time difference of arrival (TDoA) measurements. RSSI, AoA, RTT, ToA and/or TDoA measurements may be processed in any appropriate way (e.g. by filtering, averaging or otherwise) as part of determining the data representative of distances, e.g. in order to improve the accuracy or quality of the location determination process. The processing system may comprise a networked server, which may be distinct from said devices of the network (i.e. the devices whose locations are involved in the system of linear equations). However, in other embodiments, the processing system may comprise a processor of one or more of the devices—e.g. one of the mobile devices. In some embodiments, the location system (i.e. positioning system) may comprise one or more of the devices. It may comprise all of said devices of the network. (It will be appreciated, however, that these devices may in some embodiments be part of a larger network, comprising other devices that are not included in the system of linear equations—e.g. devices for which no distance data, or insufficient distance data, is received.) It may be a real-time location system (RTLS). The processing system may receive the distance data, e.g. from one or more of the devices of the network, or it may be configured to generate some or all of the distance data, e.g. by processing measurement data, such as RSSI values, received from one or more of the devices of the network. Receiving the data may comprise accessing the data within the processing system, e.g. reading the data from a memory of the processing system. The values representative of distances between pairs of devices may be noisy or approximate data. Processing the values to determine the coefficients may comprise scaling the values based on a level of noise in the set of values. The system of linear equations may be determined, at least in part, by said set of coefficients. The set of coefficients may comprise a respective coefficient (which may be a zero or non-zero value) for each pair of devices of the network of devices. Each coefficient may depend on and/or be calculated using one or more of said set of values representative of distances between respective pairs of devices of the network. In some preferred embodiments, each of the coefficients is a respective linear combination of one or more of the values of the set of values. This can ensure the system of linear equations remains self-consistent and that a solution can be found. The network of devices may comprise one or more devices having a known location— i.e. for which two-dimensional or three-dimensional coordinates are determined or received before the system of linear equations is solved. The device or devices of known location may comprise mobile devices and/or fixed devices. The known location may be exact or may be an estimate. In some embodiments, the locations of one or more of the mobile devices may be determined, or may have been determined, using other positioning means, such as from GPS measurements. In some embodiments, the network of devices comprises one or more of the devices having a fixed location, such as a static beacon configured to exchange wireless signals with the mobile devices. In some embodiments, the network of devices comprises one or more mobile devices, each configured to determine a position of the mobile device (which may be used as a “known” position) using a positioning method that does not use wireless signals transmitted between the mobile device and other devices of said network of devices—e.g. using a global navigation satellite system (GNSS). Devices of known location may be used to provide an absolute position reference and may be used to determine a respective location of each of the set of one or more mobile devices having unknown location (i.e. devices having locations that are determined by solving the system of linear equations). In particular, a respective two-dimensional or three- dimensional coordinate of one or more devices of known location may be used for determining the two-dimensional or three-dimensional coordinates representative of the respective location of each of the set of one or more mobile devices—e.g. used for determining the set of coefficients for the system of linear equations and/or used for solving the system of linear equations. In some embodiments, the method comprises determining the location of a single mobile device of the network of devices. The single mobile device may, in some embodiments, be the only mobile device of the network of devices (i.e. with the other devices having fixed locations). However in some embodiments, the method comprises determining the location of a plurality of mobile devices of the network of devices. The network of devices may comprise exclusively mobile devices, or may comprise one or more fixed devices in addition to the set of one or more mobile devices. In embodiments in which two-dimensional coordinates are determined, the two- dimensional plane may correspond to a physical environment in which the network of devices are located. It may be a horizontal plane. The two-dimensional plane may correspond to a two-dimensional view of, or area in, the real world, which may be an indoor environment, such as a room or a floor of a building. The method may comprise using the determined 2D coordinates to identify one or more respective two- or three- dimensional positions in the environment, e.g. such that a position of each of the set of one or more mobile devices within the room or on the floor of the building may be determined. In embodiments in which three-dimensional coordinates are determined, the three- dimensional space may correspond to a physical environment in which the network of devices are located. The three-dimensional space may correspond to a three- dimensional volume in the real world, which may be an indoor environment, such as a room or a floor of a building or an entire building. The method may comprise using the determined 3D coordinates to identify respective positions in the environment, e.g. such that a position of each of the set of one or more mobile devices within a room or building may be determined. One or more of the 2D or 3D coordinates, or positions corresponding to the coordinates, may be output, e.g. as data over a network, or rendered visually on a display screen, such as on a map or plan of an environment in which the network of devices are located. The location system may comprise a display screen for displaying data representative of one or more of the determined coordinates. The system of linear equations may comprise a plurality of variables corresponding to coordinates of devices of unknown location. It may additionally comprise one or more constant terms corresponding to coordinates of respective devices of known location (although this is not essential in all embodiments or situations). In some embodiments in which two-dimensional coordinates are determined, the two- dimensional coordinates may be represented as complex values in the complex plane. The two-dimensional coordinates for respective devices may be represented in the system of linear equations by complex values or complex variables. This may allow efficient algorithms to be used to solve the system of linear equations, for example by finding complex eigenvalues. The system of linear equations may be homogeneous or inhomogeneous. For example, it may be homogeneous if there are no devices of already-known location, and inhomogeneous if coordinates for one or more devices are known. The system of linear equations may be solved in any appropriate way. In some embodiments, solving the system of linear equations comprises performing singular value decomposition, or performing Gaussian elimination, or performing lower–upper (LU) decomposition. The set of coefficients may be stored as data representing a square matrix. The square matrix may comprise a respective row and a respective column for each device of the network of devices. In some embodiments, solving the system of linear equations may comprise calculating an inverse matrix. In such embodiments, the two-dimensional or three- dimensional coordinates may be determined from the product of the inverse matrix and a vector. In such embodiments, the vector may comprise one or more coordinates corresponding to the known locations of one or more devices of the network. The set of coefficients may be filtered before solving the linear system. It may be filtered by attenuating—e.g. reducing or setting to zero—coefficients, for one or more pairs of devices, that do not meet an inclusion condition. The inclusion condition may specify or determine a maximum number of non-zero coefficients, which may depend on the number of devices or the number of mobile devices (optionally in addition to depending on other factors). This may allow poorer-quality distance data between pairs of devices to be removed, thereby improving the accuracy of the determined mobile device locations. In some embodiments, at least one coefficient is set to zero. By reducing the number of non-zero terms in the system of linear equations, the efficiency with which the system of linear equations can be solved may be improved— e.g. by increasing the sparsity of a matrix representation of the system. In some embodiments in which two-dimensional coordinates are determined, determining the set of coefficients may comprise, for at least one, or each, device, identifying one or more sets of three further devices of the network. It preferably comprises identifying one or more sets of three further devices that form a respective triangle surrounding said device, in the two-dimensional plane. It may comprise using distances between said device and the three further devices, determined from the distance values, when determining the set of coefficients. Determining the coefficients may comprise performing a bilateration process (e.g. relating to the intersection of circles in two dimensions). It may comprise determining a local coordinate system determined for said device and determining coefficients that represent the locations of the three further devices in said local coordinate system. Determining the set of coefficients may comprise excluding, from the system of linear equations, distances between said device and any further device that is not part of a set of three devices that satisfies a geometric filter for said device. In some embodiments, in which two-dimensional coordinates are determined the geometric filter for a device may only pass sets of three further devices that form a triangle that surrounds the device or that is proximate to the device—e.g. for which the distance between said device and a point on the triangle that does not surround the device is within a threshold distance. The threshold distance may be determined based on an uncertainty in the determined distances between pairs of devices of the triangle. Determining the coefficients for the linear system may comprise applying a geometric filter based on Heron’s rule. In some embodiments in which two-dimensional coordinates are determined, coefficients that are associated, in the linear system, with locations of devices forming the triangles may be weighted by a variance determined from an analysis of triangles in the 2D plane. Such a variance may provide a measure of the uncertainty of the relative device locations. In some embodiments a variance is determined for a device surrounded by three further devices, based on a comparison between the area of the triangle formed by the three further devices and the sum of the areas of the three triangles formed by said device and each respective subset of two of the three further devices. In some embodiments, one or more sets of three devices forming a triangle surrounding a device may be excluded based on such a variance. In some embodiments, a plurality of sets of three further devices may be sorted by variance. In some embodiments, only distance values for sets of devices up to a threshold number of (e.g. the top four) sets of devices, ranked by variance, may be considered when determining the coefficients. In some embodiments, for each device, when identifying sets of three devices forming a triangle surrounding said device, the devices may be limited to include only up to a threshold number of (e.g. up to six) other devices. Devices may be selected by favouring devices having a known location and/or lower variance. In some embodiments in which three-dimensional coordinates are determined, determining the set of coefficients may comprise, for at least one, or each, device, identifying one or more sets of four further devices of the network. It preferably comprises identifying one or more sets of four further devices that form a respective tetrahedron surrounding said device, in the three-dimensional space. It may comprise using distances between said device and the four further devices, determined from the distance values, when determining the set of coefficients. Determining the set of coefficients may comprise performing a trilateration process (e.g. relating the intersection of spheres in three dimensions), which be done in respective of each such set of four further devices. The trilateration process may, in some embodiments, comprise a two-dimensional bilateration operation (e.g. to determine a 2D position of one of the further devices relative to another two of the further devices). It may comprise determining a local coordinate system determined for said device and determining coefficients that represent the locations of the four further devices in said local coordinate system. In some embodiments in which three-dimensional coordinates are determined, determining the set of coefficients may comprise excluding, from the system of linear equations, distances between said device and any further device that is not part of a set of four devices that satisfies a geometric filter for said device. In some such embodiments, the geometric filter for a device may only pass sets of four further devices that form a tetrahedron that surrounds the device or that is proximate to the device—e.g. for which the distance between said device and a point on the tetrahedron that does not surround the device is within a threshold distance. The threshold distance may be determined based on an uncertainty in the determined distances between pairs of devices of the tetrahedron. In some embodiments in which three-dimensional coordinates are determined, coefficients that are associated, in the linear system, with locations of devices forming the tetrahedra may be weighted by a variance determined from an analysis of tetrahedra in the three-dimensional space. Such a variance may provide a measure of the uncertainty of the relative device locations. In some embodiments a variance is determined for a device surrounded by four further devices, based on a comparison between the volume of the tetrahedron formed by the four further devices and the sum of the volumes of the three tetrahedra formed by said device and each respective subset of three of the four further devices. In some embodiments, one or more sets of four devices forming a tetrahedron surrounding a device may be excluded based on such a variance. In some embodiments, a plurality of sets of four further devices may be sorted by variance. In some such embodiments, only distance values for sets of devices up to a threshold number of (e.g. the top four) sets of devices, ranked by variance, may be considered when determining the coefficients. In some embodiments, for each device, when identifying sets of four devices forming a tetrahedron surrounding said device, the devices may be limited to include only up to a threshold number of (e.g. up to six) other devices. Devices may be selected by favouring devices having a known location and/or lower variance. More generally, in some embodiments, for each mobile device of the set of one or more mobile devices, the number of non-zero coefficients associated with the location of the mobile device in the system of linear equations may be limited by a predetermined maximum threshold number—comprising no more than six non-zero coefficients, irrespective of the number of mobile devices in the network. This can make the system more efficient to solve, and may improve accuracy if data of lower quality is excluded. In some embodiments in which two-dimensional coordinates are determined, the method may nonetheless be used to determine the location of each of the set of one or more mobile devices in three dimensions. Determining the third coordinate (e.g. corresponding to height or altitude) may be done by a different mechanism—e.g. by selecting the set of one or more mobile devices from a larger group of devices, using a signal strength condition that limits all the devices to being on the same floor of a building (which may then all be assigned a common height value, corresponding to a known height of the floor). In some embodiments, it may comprise solving a further system of linear equations, determined by the set of coefficients, to determine three- dimensional coordinates representative of a respective location of each of the set of one or more mobile devices in three dimensions. The processing system may comprise one or more processors. It may comprise a single processing device—e.g. a workstation located in the same environment as the mobile devices—or a network of processing devices—e.g. a distributed cloud server accessed over the Internet. It may comprise a memory for storing data and said software. It may comprise one or more wired or wireless network interfaces. It may comprise a display and/or other input-output peripherals. Features of any aspect or embodiment described herein may, wherever appropriate, be applied to any other aspect or embodiment described herein. Where reference is made to different embodiments or sets of embodiments, it should be understood that these are not necessarily distinct but may overlap. BRIEF DESCRIPTION OF THE DRAWINGS Certain preferred embodiments of the invention will now be described, by way of example only, with reference to the accompanying drawings, in which: FIG.1 is a schematic drawing of a peer to peer positioning system embodying the invention; FIG.2 is a schematic drawing of a static transmitter unit, a mobile device, and a server for use in the peer to peer positioning system; FIG.3 is a flow chart of operations carried out by the server to estimate the position of the mobile device in a two-dimensional plane; FIG.4 is a schematic diagram of a configuration of devices in a two- dimensional plane determined by the peer to peer positioning system; FIG.5 is a schematic diagram illustrating a process of determining relative distances between devices of the peer to peer positioning system in a local coordinate system; and FIG.6 is a schematic diagram illustrating a process of determining the position of a device of the peer to peer positioning system in terms of weighted positions of other devices of the peer to peer positioning system; FIG.7 is a schematic diagram illustrating a process of selecting subsets of devices to use for determining the position of a particular device of the peer to peer positioning system; FIG.8 is a simplified example of a system of linear equations solved by the server to estimate the positions of devices; FIG.9 is a flow chart of operations carried out by the server to estimate the position of the mobile device in a three-dimensional space; and FIG.10 is a schematic diagram of a configuration of devices in a three- dimensional space determined by the peer to peer positioning system. DETAILED DESCRIPTION FIG.1 shows elements of a peer-to-peer positioning (i.e. location-determining) system that may be used in, for example, a shopping mall, in order to determine the locations of shoppers within the shopping mall, by determining the locations of mobile devices, such as smartphones or electronic tags, carried by one or more of the shoppers. Of course, this is just one example environment, and the positioning system could also be used in warehouses, hospitals, domestic homes, vehicles, or many other indoor or outdoor environments. FIG.1 shows a top-down view of a single-storey environment 1 (e.g. a room, or a floor of a building) in which a peer to peer positioning system 100 is operating. Other details, such as people, walls, doors, furniture, etc. are omitted for simplicity. The system 100 provides a networked server 2, configured to communicate with a plurality of mobile devices 3, such as smart phones or mobile tags, and static beacons 4. The static beacons 4 are distributed throughout the environment 1 at fixed locations—e.g. fastened to walls or to the ceiling. The mobile devices 3 may move throughout the environment 1, for example being carried by people or fastened to mobile equipment. The mobile devices 3 are configured to transmit and receive wireless signals, such as Bluetooth Low Energy™, WiFi™ or ultrasound signals, or a combination of different signal types, between one another and to and/or from the static beacons 4. The mobile devices 3 and the static beacons 4 (collectively “devices”) are further configured to communicate with the server 2, which may be direct communication or which may be via other devices—e.g. operating as a wireless mesh network, having any appropriate topology. These components cooperate to provide a peer to peer positioning system 100, configured to estimate the positions of the mobile devices 3 as two-dimensional (x, y) coordinates within the single-storey environment 1 or as three-dimensional (x, y, z) coordinates. Based on data received at the server 2 relating to wireless signals exchanged by the mobile devices 3 and the static beacons 4, the server 2 is configured to determine locations of the mobile devices 3 as will be described in the following. Although in this embodiment the server 2 is shown within the same environment 1 as the devices 3, 4, in other embodiments it could be located remotely, e.g. at an Internet server farm, and be in communication with one or more of the devices 3, 4 via a network gateway or base station. While the embodiment shown in Figure 1 shows a plurality of mobile devices 3, it will be appreciated that the peer to peer positioning system may also or alternatively be used to determine the location of a single mobile device 3, without requiring signal transmissions between the mobile device 3 and any other mobile device of the system, provided that a sufficient number of static beacons 4 are located in the vicinity of the mobile device 3. Thus, the term “peer to peer”, as used herein, should be understood as encompassing “device to device” positioning between devices that may be dissimilar from each other, e.g. between a mobile device 3 and one or more static beacons 4. FIG.2 shows a representative mobile device 230 of the mobile devices 3, a representative beacon 240 of the static beacons 4, and the server 2 in more detail. The mobile device 230 has a controller 231 comprising a processor and one or more radio modules, a battery 233 for supplying power to the mobile device 230, and a memory 232 for storing data and software 10 for the controller 231. The mobile device 230 also contains an antenna 234, for transmitting and receiving radio signals between the mobile device 230 and other mobile devices 3 and beacons 4 of the peer to peer positioning system 100. The beacon 240 has a controller 241 comprising a processor and one or more radio modules, a battery 243 for supplying power to the beacon 240, and a memory 242 for storing data and software for the controller 241. The beacon 240 also contains an antenna 244, for transmitting and receiving radio signals between the beacon 240 and other beacons 4 and mobile devices 3 of the peer to peer positioning system. The beacon 240 may also have a wired interface, e.g. for connecting the beacon 240 to the server 2 over a wired Ethernet or fibre optic connection. The networked server 2 comprises a processor 227 and a memory 222 for storing data and software. It may be static and externally-powered. It may comprise an antenna 224 for direct radio communication with one or more of the mobile devices 3 and beacons 4 of the peer to peer positioning system, or it may have a wired network connection to one or more of the beacons 4, or both. The networked server 2 may be a single physical device, or may be a distributed (e.g. cloud) server system. It may be in the same building, or on the same site, as the mobile devices 3 and static beacons 4, or it may be remote from these components, e.g. in a remote data centre accessed over the internet. While the devices shown in FIG.2 comprise respective antennae, they may instead transmit and receive wireless signals other than radio signals, such as ultrasound signals, using acoustic transducers (i.e. loudspeakers and microphones). In use, each mobile device 3 of the peer to peer positioning system 100 is configured to listen for wireless signals from other mobile devices 3 and beacons 4 of the system 100, and to send wireless messages at regular or irregular intervals, using its antenna 234. Such communications can be used to estimate inter-device distances between pairs of the devices 3, 4 that are in effective communication range of each other, using conventional techniques, such as from received signal strength indicators and/or time- of-flight measurements. For example, distance estimates from a mobile device 3 to one or more other devices 3, 4 may be determined based on one or more received signal strength indicators (RSSI) of BLE advertising messages received by the device 3. In other examples, they could be determined through Bluetooth angle of arrival (AoA), WiFi RSS or round trip time (RTT), ultrawideband (UWB) RTT, ultrasound time of arrival (ToA), ultrasound time difference of arrival (TDoA) measurements, or in any other appropriate way. Such measurements may be processed in any appropriate way in order to determine these inter-device distance estimates (e.g. filtered, averaged, adjusted, combined, etc.). Inter-device distances between the mobile devices 3 and with the static beacons 4 are determined and stored in the memories 232 of the mobile devices and/or the memories 242 of the static beacons as inter-device distance data. The inter-device distance data is transmitted to the server 2, where it is saved in the memory 222 of the server, and processed by the processor 227 to determine locations of the mobile devices 3. Inter-device distance data may be collected by one or more of the devices 3, 4 and transmitted to the server 2 at intervals, or may be transmitted each time an interaction is recorded between mobile devices 3 and/or static beacons 4. A database of the inter-device distance values between the mobile devices 3 and/or the static beacons 4 is therefore built up and stored in the memory 222 of the server 2, and is updated dynamically as new inter-device distance measurements are recorded. As the positions of the static beacons 4 are fixed, their absolute locations can be known to the server 2, as (x, y) or (x, y, z) coordinates relative to a predetermined origin and orientation, and hence do not need to be determined using the position determination process. This information may be transmitted to the server 2 during transmission of the inter-device distances, or may be stored in the database at the time the system is set up. The absolute position of certain mobile devices 3 may also be known, for example based on GPS measurements obtained by the devices 3. For example, considering a mobile device 3 initially located outside of the indoor environment shown in Figure 1, the position of the mobile device 3 may be known based on GPS measurements. When the mobile device 3 enters the indoor environment shown in Figure 1, and GPS signal is lost, it may begin operating as part of the peer to peer positioning system 100, and may transmit its last position before the GPS signal was lost to the server 2. This may then serve as an absolute position reference for the mobile device 3—i.e. providing a “known” position. Similarly, if the position of a mobile device 3 is known from other sources, its known position may be transmitted to the server 2, e.g. during transmission of the inter-device distances. Known positions serve as reference locations from which the absolute locations of mobile devices 3 of unknown position can be determined. Thus, in some embodiments, the peer to peer positioning system 100 may be able to determine the locations of devices of unknown position based only on the known positions of one or more mobile devices 3, i.e. without necessarily requiring the provision of any static beacons 4 of fixed position. Having established a database of inter-device distances, the peer to peer positioning system 100 can process the distance values to determine a constellation of the positions of some or all the mobile devices 3 of the system 100 in two- or three- dimensions as will be described in the following. The approach described herein allows the position of mobile devices 3 to be determined indirectly, based at least in part on position estimates calculated for other mobile devices 3. Thus even mobile devices 3 that are not within wireless range of any static beacon 4 can potentially be accurately positioned. An overview of the position determination process according to a first set of embodiments is illustrated in Figure 3, which shows a flow chart of operations carried out by the server 2 to estimate the positions of some or all of the mobile devices 3 of the system 100 in a two-dimensional plane. In step 301, the server 2 identifies all the mobile devices of unknown position. The server 2 then performs a looped set of steps for each of these devices of unknown location in turn. In step 303, the server 2 determines, from the inter-device distance values, all combinations of three other devices j (which may be mobile devices 3 or static beacons 4) for which there is recent distance data (i.e. data recorded within a threshold time period, such as 15 seconds) between the device i and the other device j. In step 305, a reduced subset of combinations of three devices j is determined, which are used in subsequent steps to determine the position of the mobile device i. This subset is chosen based on a condition that the mobile device i should be located within or sufficiently close to a triangle formed by the combination of devices j, in the (x, y) plane. This determination is made using a criteria based on Heron’s formula as will be explained below. In step 307, a bilateration process is used by the server 2 to determine, from the inter- device distance data, the relative positions of the mobile device i and each of the surrounding devices j contained in the subset of combinations, in a local coordinate system. In step 309, the position z _i of the mobile device i is determined in the form of a linear weighted sum of the positions z _j of the devices j, i.e. in the form z _i =w _ij z _j , for real-valued coefficients w _ij , with the (x,y) positions represented as points in the complex plane by complex values z _i and z _j . These steps 303 to 309 are repeated for each device i of unknown location. Then, in step 311, data (e.g. a matrix of coefficients) representing a system of linear equations is generated, in which the position of every device 3, 4 in the system 100 is expressed in terms of a respective weighted sum of positions of relevant other devices. The matrix may optionally be filtered to increase its sparsity. In step 313, the terms of the matrix relating to devices 3, 4 with already-known positions are multiplied out, in order to reduce the size of the matrix by the number of known device positions. In step 315, the system of linear equations is solved (e.g. by inverting the remaining matrix), to determine a vector containing the positions of the mobile devices 3 of previously unknown position in a two dimensional plane. Having outlined the process by which the position determination process is performed by the server 2, the steps of the process according to this first set of embodiments will now be described in more detail with reference to FIGs 3-7. Within the first step 301, all devices (i.e. mobile devices 3 and beacons 4) to be positioned within an area or environment are identified. This may be all of the devices located on a floor, or in a room of a building, or may be all of the devices of the peer to peer positioning system 100. The selection of devices within a room or limited to a floor of a building may be made based on a level of attenuation of signals transmitted between devices. High-frequency signals in the ISM band (2, 4 and 5 GHz) experience very high attenuation when crossing walls and floors, and hence devices in a common room or floor can be easily grouped. Thus a third-dimension of height (e.g. which storey a device is on) may be determined before the present two-dimensional positioning process is initiated, such that two-dimensional (x, y) coordinates in the horizontal plane may also be sufficient for resolving the location of a device in three dimensions (x, y, z) if required. Alternatively, as will be described below, a second set of embodiments use a variant process that can directly determine device positions in three-dimensions, from the wireless signals. An inter-device distance graph for all of the selected devices is already determined by data stored in a database of the server 2, including the inter-device distances of the selected devices, and an indication of whether the position of each of the devices is known or unknown. As described above, known device positions include the position of the static beacons 4 and the position of any mobile device 3 whose absolute position is known from other means, for example from GPS measurements. These locations may be determined by latitude & longitude, or relative to a predetermined local coordinate system for the environment 1, or in any other appropriate way. Using the inter-distance device graph, the relative positions of the selected devices are then determined as will be described in the following. The positioning method according to this first set of embodiments, as described below, exploits the fact that, if the devices are considered to be located on a 2D plane (which is a reasonable assumption in many cases, such as when the devices are known to all be on the same floor of a building), the position of each device can be expressed as a linear combination of the 2D positions of the devices that surround it. The resulting system of linear equations can be solved simultaneously, to determine estimates for a plurality of devices simultaneously. The surrounding devices may all be mobile devices or may all be fixed beacons, or may be any combination of mobile and fixed devices. In step 303, for a given mobile device i with unknown position z _i , the server 2 determines every combination of three other devices j. This may be achieved using conventional methods for finding all C(N, 3) (i.e. “N choose 3”) sets, such as using a Cartesian product list. In step 305, these are filtered to find only those combinations having positions z _j forming a triangle surrounding, or nearly surrounding (i.e. close to), the position z _i . This is done by applying a geometric filter based on Heron’s rule. An exemplary combination of three devices j identified in this way is shown in FIG.4. FIG.4 shows a mobile device i at a position z _i (represented by coordinates (x,y) in the 2D plane) surrounded by three mobile devices j ₁ , j ₂ , j ₃ located at positions z _j1 , z _j2 , z _j3 respectively (and represented by coordinates (x ₁ ,y ₁ ), (x ₂ ,y ₂ ), (x ₃ ,y ₃ ) in the 2D plane). The six pairwise inter-device distances between each of the four mobile devices, determined by wireless ranging and stored in the memory 222 of the server 2, are shown in FIG.4 as distances a-f. These inter-device distances a-f are used to determine whether the mobile device i lies within, or close to, a triangle defined by the positions zj of the devices j1, j2, j3 in the x-y plane, by considering the respective areas of triangles formed by the device positions. The area of the largest triangle formed by the devices j1, j2, j3, shown in FIG.4 with edges a-b-c, and the areas of the three inner triangles, having edges a-d-e, d-c-f and e-f-b respectively, can be calculated using Heron’s rule, which states that the area A of a triangle having sides of lengths p, q and r, equals: where s is the semi-perimeter of the triangle, given by: Using the above equations, the server 2 determines whether the three devices j ₁ , j _2, j ₃ of every possible combination surround the mobile device i by comparing the sum of the areas of the inner triangles a-d-e, d-c-f, and e-f-b (shown as A1, A2, A3 in FIG.4) to the area A of the outer triangle a-b-c. If the device i is contained within the triangle, these two values will be equal, or approximately equal when allowing for measurement errors. As the inter-device distances determined using RSSI measurements in the peer to peer positioning system 100 have some associated error, a relaxed comparison of the areas is made according to: _ where is selected based on properties of the peer to peer positioning system, such as an expected RSSI error, and is typically in the range 1 ≤ By applying this geometric condition to the distance data for each combination of three devices, all the combinations that effectively surround the mobile device i can be determined. For each combination of devices that satisfy the above condition, a measure of uncertainty of the relative device positions may be determined by comparing the variance of the inner triangle areas to the mean area The results of this comparison may be used to determine a reduced subset of combinations of three devices to be used in later calculations, e.g. by excluding combinations of devices which form triangles for which the uncertainty is high, such as above a threshold. In some but not all embodiments, it may be used by the server 2 to influence the calculation of the coefficients w _ij in the linear weighted sum. Having determined all possible combinations of three devices j ₁ , j _2, j ₃ surrounding the mobile device i for which a complete inter-device distance set is available (i.e. having six known inter-device distances between the four devices involved), the server 2 then, in step 307, determines the relative positions of the devices i and j1, j2, j3 in a local coordinate system. The position of the device i is determined within a local reference frame defined by the devices j1, j2, j3 by performing a series of bilaterations—that is, by selecting an arbitrary reference frame in which the position of one of the devices, e.g. j2 in FIG.4, is equal to (0,0) and the position of another of the devices, e.g. j3 in FIG.4, is equal to (b,0), the position of the final device, e.g. j1, can be calculated in this local coordinate system from the known distances a, b, c using bilateration (i.e. based on the intersection of circles in the plane). Further bilaterations allow the position of the mobile device i to be similarly determined. For a triangle of arbitrary coordinates (x1,y1), (x2,y2), (x3, y3), with distance r1 between (x1,y1) and (x2,y2), and distance r2 between (x1,y1) and (x3,y3), the position of (x1,y1) is given by where FIG.5 shows an example of applying this to the outer triangle from FIG.4, to express the position of the device j1 as (x1,y1)=(l, h), where r1 is the measured distance between j ₂ & j ₃ (equal to “a” from FIG.4), r ₂ is the distance between j ₁ & j ₃ (equal to “c” from FIG.4), and D is the distance between j2 & j3 (equal to “b” from FIG.4). Once the positions of the devices j1, j2, j3 are determined in a local coordinate system, the bilateration process is repeated for each of the three inner triangles to determine respective estimates of the position of the mobile device i, located within the triangle formed by devices j ₁ , j _2, j ₃ in the local coordinate system. While successive bilaterations could lead to increased complexity, as two solutions (+h,-h) are possible, by ensuring that a consistent order of the devices is maintained when assessing subsequent bilaterations, the positive h solution can be consistently selected. The three bilateration results are rotated using the orientations of the respective baseline sides to align with the original local reference frame defined by the edges a,b,c. The three bilateration position estimates obtained for each sub-triangle will typically differ due to measurement noise. This provides a good way of assessing the variance of the position (pos_variance) estimate for the mobile device i by calculating the second norm with respect to the mean position. This true positional variance may be used, in some embodiments, to weight the coefficients by a factor 1/pos_variance. Having determined a position estimate z _i for the mobile device i in the local coordinate system defined by the triangle a-b-c, its position z _i can be mapped back to the positions zj of the devices j1, j2, j3 by the coefficient factors wij. The position zi of the mobile device i can be expressed in terms of the device positions j1, j2, j3 in vector notation. The theory behind this is illustrated with reference to FIG.6. FIG.6 shows a triangle formed from devices j1, j2, j3, within which a mobile device i is located. The position P of the mobile device i can be expressed as a two-dimensional vector in terms of the positions P1, P2 and P3 (equivalent to zj1, zj2 and zj3 respectively) of the devices j1, j2, j3. The server 2 may, in some embodiments, be configured to perform calculations on these positions by representing them as complex numbers in the complex plane. Specifically, the position P of the mobile device i can be written in terms of a vector v from v0 to the position P of the mobile device i as: where v ₁ and v ₂ are the vectors from v ₀ to the two other vertices of the triangle, and where α and β are constants. Solving for α and β gives: where d is the determinant of the matrix formed by the column vectors u and v given by: By calculating ^ and ^, and writing the vectors v1 and v2 in terms of the vertex positions P1, P2 and P3 as v1=P1-P2 and v2=P3-P2, the position P of the mobile device i can be determined as Taking the real-valued coefficient wij for each device j1, j2, j3 can be calculated by the server 2 in terms of the constants ^, ^ and ^, such that the position of the mobile device i can be expressed in terms of the positions of the surrounding devices j1, j2, j3. The above calculation of also provides another means for checking that the position of the mobile device i is inside of the triangle formed by the devices j, which is true only In some embodiments, if this relationship is not satisfied, the results are discarded. Following the above calculations, ^, ^ and ^ are the normalised weights for each triangle's node positions. These can be further weighted using the positional variance found for the specific triangle, as described above, such that the coefficients are further weighted by a factor of _ The coefficients are therefore able to account for positioning inaccuracy of the devices. Further weighting of the coefficients may be applied based on the recency (i.e. freshness) of the distance data used to determine the positional variance. This may be achieved by multiplying the positional variance by a pre-factor that decreases linearly with time up to a threshold time (e.g. 15 seconds), after which the pre-factor is reduced to zero. Additional sources of positioning error may optionally be taken into account by further modifying the weighting of the coefficients. For example, in the case where the position of one or more of the devices is “known”, e.g. from a GPS measurement, its “known” position may still have an associated error due to positional measurement inaccuracy. In this case, the uncertainty in the device’s position may be taken into account, e.g. by adding a variance term, based on the positional measurement inaccuracy, to the positional variance To estimate the impact of positional variance from one or more of the devices j on the variance of device i, the variances may be propagated using the method of partial derivatives for each triangle, such that: where F is the functional relationship between zi and zj. The above processes allow the server 2 to determine the position z _i of a mobile device i as a consistent linear combination of the positions z _j of devices j that surround it, each weighted by a coefficient w _ij . The applicant has recognised that by determining coefficients using the process described above, a self-consistent system of linear equations can be established that relates the positions of devices of the peer to peer positioning system to each other. If the 2D coordinates of the devices are considered as complex numbers in the complex plane, and represented as a column vector for the set of all M devices (including those of known and unknown positions), the system can be expressed as a 2D matrix equation, with matrix elements corresponding to the coefficients wij, such that If the diagonal terms are set such that then this is equivalent to which can be expressed in matrix form as If the server 2 can determine sufficiently many non-zero coefficients for pairs of devices forming part of the peer to peer positioning system 100, the position of every device zj can be simultaneously calculated through conventional techniques for solving linear systems. The server 2 may be configured to solve the system of linear equations by any appropriate method. If the position of some devices 3, 4 is known, the terms of the matrix relating to any devices with known positions, zk, may be multiplied out by the server 2 prior to solving the linear system, in order to reduce the size of the matrix by the number of known device positions. For example, in a system comprising N mobile devices of unknown position and k devices of known position (e.g. static beacons and/or mobile devices having known positions), the matrix equation may be represented as Because the device positions zN+1 to zN+k are known, these can be separated out as The known positions zN+1 and zN+k and the calculated coefficients relating to the devices of known position can be expressed as a constant vector b, leading to an inhomogeneous matrix equation, where the values of b are calculated according to The server 2 can find a solution to this inhomogeneous matrix equation, for the unknown positions z1 to zN, using conventional techniques, such as singular value decomposition, Gaussian elimination, or lower–upper (LU) decomposition to determine the unknown positions of the devices z ₁ to z _N . FIG.8 provides a simplified example of such a system of linear equations, showing a matrix, w _ij , of real-valued weight coefficients, a vector z _j of unknown device positions (where each z _j is the complex number formed from the device’s Cartesian coordinates as x _j + i.y _j ), and a vector b _i which is calculated by the server 2 as -sum(w _ik z _k ) from the positions, z _k , of k known devices. The server 2 may combine as many triangles as possible to populate the matrix w _ij . However, for high device counts it is advantageous to keep the matrix wij sparse, such that less computationally intensive methods for solving the linear system (e.g. for performing matrix inversion) may be used. It may therefore be beneficial, at least in some embodiments, for the server 2 to filter the matrix by setting some coefficients to zero (or otherwise attenuating them). It may do this by imposing a maximum number of triangles to be used for determining the coefficients for each device i, such as four triangles. It may additionally or alternatively impose a minimum number of non-zero matrix elements (coefficients wij), such as a minimum of four—i.e. it may keep adding triangles until the minimum is reached for each row. It may additionally or alternatively impose a maximum number of non-zero matrix elements, such as a maximum of six, that are allowed in each row. By applying both a minimum and maximum number of allowed non-zero elements per row, computational efficiency can be optimised while ensuring a sufficient number of device counts to be able to effectively determine device positions. An example of such a filtering can be explained with reference FIG.7, which shows a mobile device i of the peer to peer positioning system 100 surrounded by seven other devices j including four static beacons (shown as triangles in FIG.7) and three other mobile devices (shown as empty circles in FIG.7). In FIG.7, three triangles 701-703 have been selected by the server 2 to be used to determine the coefficients for the system of linear equations in respect of the mobile device i. FIG.7 also illustrates how the Heron criterion may, in some situations, allow a triangle 701 to be used even when the mobile device i actually lies just outside the triangle 701, because of the error margin. When selecting a subset of triangles to be used, two classes of triangles can be set: a first class of those triangles having at least one apex defined by a beacon 4 or device 3 of known location, and a second class of triangles which do not have this. Thus, triangles may be formed exclusively from beacons 4 or exclusively from mobile devices 3, or by any combination of beacons 4 and mobile devices 3. Preferably, when selecting a reduced subset of triangles to be used (e.g. restricting to four triangles), the server 2 favours selecting triangles from the first class. For devices in the same class, it may be configured to select triangles in order of lowest positional variance. This typically favours the surrounding triangles with the smallest sides. Devices 3, 4 that do not have any valid surrounding triangle of nodes or have no nearby triangles to support extrapolating are removed from the coefficient matrix wij and do not have their positions determined with the method described above. However, in some cases where a mobile device i is not surrounded by any triangle of other devices, the server 2 may be able to include it by a process of extrapolation. It may attempt to add the mobile device to the linear system of device positions using a trilateration result from the closest triangles. The coefficients wij can be constructed in the same way as when a point is inside a triangle, but extrapolated rather than interpolated. One complexity that does arise when using extrapolation is that, in contrast to when a node is inside a triangle, each of the three bilaterations from a nearby triangle results in two valid positions, resulting in a total of 2 ³ =8 different combinations. To select the true estimated position, the position error for each combination may be calculated and the position with the smallest associated error selected. While the method described in relation to Figures 3 to 8 is sufficient for the determination of device positions in two dimensions, as may be useful for determining the locations of customers within a particular floor (storey) of a shopping mall for example, the inventors have also recognised that wireless signals as disclosed herein can also be used for directly determining device positions in three dimensions. Thus a second set of embodiments will now be described, with reference to Figures 9 and 10, that performs direct 3D positioning. Similarly to the two-dimensional case described above, determination of device positions in three dimensions can be achieved using the peer to peer positioning system 100 shown in Figure 1, based on the determination of inter-device distances between pairs of mobile devices 3 and static beacons 4 and the subsequent processing of these distances using a server 2. In the three-dimensional case, inter- device distances are determined in a manner equivalent to that described above (e.g. using RSSI measurements), and are saved to the memory 222 of the server 2 in a database in the form of inter-device distance data. The positions of the static beacons may be known to the server as three-dimensional (x, y, z) coordinates relative to a predetermined origin and orientation. The peer to peer positioning system 100 can process the inter-device distance data to determine a constellation of the positions of some or all of the mobile devices 3 of the system 100 in three-dimensional space as will be described in the following. The position of each device can then be expressed as a linear combination of the 3D positions of the devices that surround it, and the resulting system of linear equations can be solved simultaneously, to determine estimates for a plurality of devices simultaneously. An overview of the three-dimensional position determination process is illustrated in Figure 9, which shows a flow chart of operations carried out by the server 2 to estimate the positions of some or all of the mobile devices 3 of the system 100 in a three dimensional space. This process is described below, with further reference to Figure 10. In the first step 901, all of the devices (i.e. mobile devices 3 and beacons 4) to be positioned within an area or environment are identified, and an inter-device distance graph for all of the selected devices (including devices having unknown positions and those whose positions are known, e.g. from GPS measurements) is generated, in a manner equivalent to that described in relation to step 301 of Figure 3. In step 903, for a given mobile device i with unknown position q _i the server 2 determines every combination of four other devices j (which may be mobile devices 3 or static beacons 4) for which there is recent distance data (i.e. data recorded within a threshold time period, such as 15 seconds) between the mobile device i and the other devices j. As described above in relation to the 2D case, this may be achieved using conventional methods for finding all C(N, 4) (i.e. “N choose 4”) sets, such as using Cartesian product list. These are filtered, in step 905, to generate a reduced subset of combinations of four devices j having positions q _j forming a tetrahedron surrounding, or nearly surrounding, the position q _i in 3D space. This is conceptually similar to the determination of surrounding triangles of the 2D case, but relies on the use of tetrahedra in place of triangles. This determination is made using a criteria based on a formula similar in principle to that of Heron’s rule as will be explained below in relation to Figure 10. Figure 10 shows a mobile device i at a position qi (represented by coordinates (x0, y0, z0) surrounded by four mobile devices j1, j2, j3, j4 located at positions qj1, qj2, qj3 and qj4 respectively (and represented by coordinates (x ₁ ,y ₁ ,z ₁ ), (x ₂ ,y ₂ ,z ₂ ), (x ₃ ,y ₃ ,z ₃ ) and (x4,y4,z4)). The mobile devices j1-j4 sit at the corner points of a tetrahedron surrounding the mobile device i in 3D space. The inter-device distances between each of the mobile devices j1-j4 are shown in Figure 10 as a1-f1, while the distances between the mobile device i and the surrounding devices j1-j4 are shown as α1, β1, δ1 and γ1. While Figure 10 shows a mobile device i surrounded by four further mobile devices j1-j4 it will be appreciated that any or all of the surrounding mobile devices j ₁ -j ₄ could instead be respective fixed beacons 4. The inter-device distances are used to determine whether the mobile device i lies within, or close to, a tetrahedron defined by the positions qj of the devices j1-j4. This is achieved by comparing the combined volume, VT, of the sub-tetrahedra formed by three of the four devices j1-j4 and the mobile device i, and having volumes V1-V4 respectively. In other words, a determination is made as to whether: In terms of the device positions, this can be expressed as where V(j1,j2,j3,j4) is the volume of the tetrahedron with vertices j1, j2, j3, j4. When comparing the combined volume V _T to the volume of the sub-tetrahedra using this equation, the equality will hold for any position inside the tetrahedron, while the combined volume of the sub-tetrahedra will exceed V _T for any position outside the combined tetrahedron. Thus, if the mobile device i is contained within the combined tetrahedron of volume VT, the two sides of the equation above will be equal, or approximately equal when allowing for measurement errors. The volume of a tetrahedron can be computed using the length of its sides, using existing formulae similar in principle to that of Heron’s rule described above in relation to the two-dimensional case. In general, for a tetrahedron with vertices p _a , p _b , p _c , p _d , the volume V _tetr of the tetrahedron can be given according to the following: Thus, by inputting the inter-device distances a1-f1 between the devices j1-j4 and the device i into the equations above, the volume VT of the tetrahedron formed by the devices j1-j4, and the volumes of each of the sub-tetrahedra V1-V4 formed by the device i and three of the devices j1-j4, can be determined and compared by the server 2. This method therefore allows determination of whether the device i is contained within the tetrahedron formed by the devices j1-j4 based only on the inter-device distances a1-f1. By applying a geometric condition to the distance data in this way for each combination of four devices, all the combinations of devices j that effectively surround the mobile device i can be determined in a similar manner to described above in relation to the 2D case. A volume factor may be applied, similarly to the area factor of the 2D case, to account for errors arising from the fact that the inter-device distances are determined using RSSI measurements with an associated error. This may allow for a more relaxed comparison to be made between the volumes than the strict equality defined above. The volume factor may therefore be selected based on properties of the peer to peer positioning system, such as an expected RSSI error. Similarly to the two-dimensional case, a measure of uncertainty in device position may be determined—in particular, the server 2 may, for each combination of surrounding devices, evaluate a measure of uncertainty in the relative device position resulting from that combination by comparing the variance of the inner sub-tetrahedral volumes to the mean volume of the tetrahedron. The server 2 and may use these uncertainty measures to determine a reduced subset of combinations of devices excluding those combinations with high uncertainty (e.g. above a threshold level). Having determined all possible combinations of four devices j ₁ , j ₂ , j ₃ , j ₄ surrounding the mobile device i for which a complete inter-device distance set is available (i.e. having ten known inter-device distances between the five devices involved), the server 2, in step 907, determines the relative positions of the mobile devices i and the devices j ₁ , j ₂ , j ₃ , j ₄ in a local coordinate system. Using the labelling system shown in Figure 10, to achieve this, the position of a first one of the devices j ₁ , j ₂ , j ₃ , j ₄ , for example j ₁ , is selected as an origin, i.e. having coordinates (x ₁ ,y ₁ ,z ₁ )=(0,0,0). The position q _j2 of a second one of the devices, e.g. j ₂ , is then defined as lying on the x-axis, such that it has coordinates (x ₂ ,y ₂ ,z ₂ )=(b ₁ , 0, 0). The position qj3 of a third one of the devices, e.g., j3, is used to define an x-y plane such that j3 lies in this x-y plane. The coordinates (x3,y3,z3) of this point are found using circle-based bilateration, to define the coordinates of the third of the devices as where ∅ is the angle between a1 and b1, and can be given by: The position of the fourth one of the devices, e.g. j4, is then defined as lying above the x-y plane with a positive height in the z-direction. The coordinates of this point may be found using sphere-based trilateration. In this process, the location of each of the devices j _1, j ₂ and j ₃ , now defined as lying in the plane z=0, is set as the centre of a sphere, the radii of which may be defined as r _j1 , r _j2 and r _j3 respectively. Two intersection points (one above and one below the z-axis) can be determined for the three spheres by defining the radii of the spheres in terms of the x, y, z coordinates of the position of the fourth of the devices, and solving for x, y and z. Although two intersection points are determined, one of these points is selected such that it satisfies the requirement that the position of the fourth one of the devices is above the x-y plane. The radius of each of the spheres can be defined using the following: which can be solved by an algorithm on the server 2 to obtain the x, y and z coordinates (x4,y4,z4) of the fourth one of the devices, j4, as Having defined the positions of the devices j1-j4 in a local coordinate system, the position of the mobile device i is then computed in this local coordinate system. This is achieved using a further sphere-based trilateration process. Specifically, a series of sphere-based trilateration operations are performed, each of which uses a different group of three of the four tetrahedron vertices as its basis. This process provides four different approximations of the position qi of the mobile device i, from which the average may be taken to more accurately determine the mobile device position. Once the positions q _j1-4 of the four devices j ₁ -j ₄ and the position q _i of the mobile device i have been estimated, vertex analysis can be performed to find the relationship between the device positions, mapping the position of the mobile device i back to the positions of the other devices j ₁ -j ₄ . The position q _i of the mobile device i is therefore determined in the form of a linear weighted sum of the positions q _j1-4 of the devices j _1-4 in step 909. This is achieved by setting up the following equation: This three dimensional linear equation can be solved to determine the coefficients b, c and d. A fourth coefficient can be calculated as: to provide a weighted connection between the point q _i and the points q _j1 -q _j4 , where w _i1 =a, w _i2 =b, w _i3 =c, w _i4 =d. The above calculation therefore allows the server 2 to determine the position q _i of a mobile device i as a consistent linear combination of the positions q _j of the four devices j that form a tetrahedron surrounding it, each weighted by a weighting coefficient w _ij defined by the coefficients a-d above. As in the two-dimensional case, the steps 903 to 909 are repeated for each device i of unknown location and set of four devices j, allowing sets of four elements with which to define a matrix W of weighting elements wij to be sequentially determined. With sufficient repletion of the tetrahedron analysis method, a system of linear equations can be generated, in which the position qi of each device i in the system is expressed as a consistent linear combination of the positions qj of devices j that surround it, each weighted by a coefficient wij. As in the two-dimensional case, this system of linear equations can be filtered, by setting some coefficients to zero or otherwise attenuating them, to increase its sparsity. This could be done by adapting any or all of the methods described above in the 2D embodiments, such as imposing a maximum number of tetrahedra to be used for determining the coefficients for each device i, such as four tetrahedra, or a minimum number of non-zero matrix elements (coefficients wij), such as a minimum of four, or a maximum number of non-zero matrix elements, such as a maximum of six, that are allowed in each row. While in the two-dimensional case described above it is convenient to express 2D positions as complex numbers, this is not as convenient in the three-dimensional case. Instead, each dimension is treated separately, resulting in three different matrix equations representing the x, y and z dimensions respectively. This is possible as the method described above uses only linear equations, such that three separate matrix equations , can be defined in step 911. The form of each of these matrix equations is similar to that described above in relation to the 2D case. For example, in a system comprising N mobile devices, the matrix equation of device positions in the x-dimension may be represented as If it is assumed, as in the two-dimensional case, that there are N devices of unknown position and K devices of known position (e.g. static beacons and/or mobile devices having known positions), the matrix equation above can be split and decomposed as described in relation to the 2D case in step 913, into an inhomogeneous matrix equation to obtain is a constant vector. This process can be applied in the same way for the y- and z-dimensions, in order to obtain the following three inhomogeneous matrix equations for the positions of unknown devices, 1 through N, in a system with K devices of known position.

Solving these three matrix equations for provides the positions of the unknown devices. Thus, in step 915, the server 2 solves the system of linear equations, e.g. by inverting the matrix or by any other applicable method, to determine a vector containing the positions of the mobile devices 3 of previously unknown position in three dimensions. As in the two-dimensional case, this operation advantageously only has to be carried out once in order to determine the positions of all of the devices of the system having unknown position. It will be appreciated by those skilled in the art that the invention has been illustrated by describing one or more specific embodiments thereof, but is not limited to these embodiments; many variations and modifications are possible, within the scope of the accompanying claims.