Leveraging Chaos for Wave-Based Analog Computation: Demonstration with Indoor Wireless Communication Signals

In sight of fundamental thermal limits on further substantial performance improvements of modern digital computational processing units, wave-based analog computation is becoming an enticing alternative. A wave, as it propagates through a carefully tailored medium, performs the desired computational operation. Yet, the necessary designs are so intricate that experimental demonstrations will necessitate further technological advances. Here, we show that, counterintuitively, the carefully tailored medium can be replaced with a random medium, subject to an appropriate shaping of the incident wave front. Using tunable metasurface reflect-arrays, we demonstrate our concept experimentally in a chaotic microwave cavity. We conclude that off-the-shelf wireless communication infrastructure in combination with a simple reflect-array suffices to perform analog computation with Wi-Fi waves reverberating in a room.

The first steps in the history of computation were analogue mechanical devices designed for specific tasks such as the Thomson brother's Harmonic Analyser to determine an input's Fourier representation [1]. With the dawn of the digital age, electronic general-purpose processors took over for good. Nowadays, the unceasing demand for ever-increasing computational power is confronted with fundamental thermal limits: already today, only a certain amount of transistors on a processor can be powered at any given time [2]. This so-called dark silicon problem recently sparked renewed interest in analogue specific-purpose computation, this time wave-based rather than mechanical, to serve as hardware accelerator. Traditional optical components are known to perform certain operations, for instance a lens yields the Fourier transform of an impinging wavefront [3,4], but remain very bulky. A recent breakthrough was the numerical demonstration that a suitably designed metamaterial block could miniaturize wave-based computing systems [5,6]. Yet, engineering such a complex structure of meta-atoms in practice with sub-micrometer tolerances remains a challenge to be mastered [7]. Consequently, recent experimental demonstrations of optical analogue computation resorted to leveraging specific effects for specific operations, such as spatial differentiation using interference phenomena arising from surface plasmon excitations [8], but do not offer a route to implementing any desired arbitrary computation operation. In parallel, the field of coherent nanophotonic circuits emerges that essentially employs a cascade of Mach-Zehnder interferometers [9,10]; despite some promising first experimental demonstrations, these devices remain inherently challenging to fabricate and scale up.
Here, we introduce the idea that in fact, rather than carefully fabricating a metamaterial or cascading interferometers, any disordered medium can be used to perform analogue computation on impinging waves as they propagate through the medium, and moreover that it may constitute a reconfigurable computation unit -subject to appropriate wavefront shaping. Our scheme disposes of any material fabrication challenges and applies to all forms of wave phenomena. To demonstrate the ease of practical implementation, we report an original proof-of-concept experiment under conditions that mimic faithfully a Wi-Fi wave in an indoor environment, and we use simple home-made phase-binary reflect-arrays to shape it. With only 96 phase-binary elements, we demonstrate a 4 × 4 complex-valued operation, as well as a 16 × 16 operation using the same physical system in a time-sequential manner.
At first sight, the complete scrambling of a wavefront as it propagates through a complex medium may seem contrarious to our objective of performing computation for which a very specific operation is to be carried out on the impinging wavefront [11]. Indeed, for a long time the random secondary sources (scatterers or reflectors) constituting a complex medium were considered undesirable and detrimental. Various novel techniques, notably Time Reversal and Wavefront Shaping (WFS), broke that paradigm by leveraging the secondary sources as degrees of freedom [12,13]. Inspired by aberration corrections with deformable mirrors in Astronomy, WFS discretizes the impinging wavefront and applies the appropriate phase and/or amplitude modulation to each segment such that after propagation through the complex medium the desired wavefront is obtained [13][14][15][16].
Initially conceived to counteract the scrambling of a wavefront in space and/or time [13,[17][18][19][20], WFS later-on paved the path towards harnessing a medium's complexity, as in the case of focusing beyond the Rayleigh limit achievable in homogeneous media [21][22][23]. Notably, WFS in complex media enabled the demonstration of programmable beam-splitters, the investigation of complex quantum-walks as well as custom-tailored mode sorting [24][25][26]. Moreover, random projections through multiply scattering media were recently shown to be a potentially powerful data preprocessor, for instance to approximate kernels in machine learning [27].
Very generally, a computation operation may be expressed as a matrix G that is applied to an input vector X, yielding the output vector Y = GX. To illustrate our general protocol, in Fig. 1 we consider the case of a 4 × 4 operation. We divide the incident wavefront into four segments, A to D, that will take the role of each of the four entries of the input vector. Moreover, we observe the wave field after propagation through the complex medium at four independent points, color-coded in Fig. 1, that will take the role of the output vector's entries. Without any shaping, each wavefront segment will yield different random contributions to each of the four observation points, as illustrated in Fig. 1(a) for the segment labeled D. More generally, if segment D corresponds to several wavefront shaper elements as in our illustration, each element has its own random contribution; stringing together all of the elements' contributions yields a random walk in the complex Fresnel representation.
By shaping the wavefront segment, these random walks can be converted into tamed walks whose resultants mimic the corresponding entries of G, as shown in Fig. 1(b) for segment D. The shaped wavefront for segment D has to simultaneously yield resultants with specific phases and amplitudes at all four observation points, a constraint that is much harder to satisfy than simple focusing. Once a library of optimised configurations for each wavefront segment is identified, an impinging wavefront can be shaped such that its subsequent propagation through the complex medium carries out the desired computational task, at the speed of wave propagation. This concept shifts the burden from designing and fabricating the physical layer (the medium) to identifying and imposing appropriate wavefronts. This is substantially simpler both in terms of the optimization problem as well as the experimental implementation. Moreover, the concept offers an open-loop reconfiguration possibility if the impact of the wavefront shaper's elements on the transmissions through the complex medium has been characterised [17,28].
Instead of encoding the input vector X into the impinging wavefront, alternatively X may be encoded directly into the wavefront shaper configuration such that a plane impinging wavefront can be used. This idea is visualized in Fig. 1(c) for X = [0 0 1 1], where segment D is in a configuration D 1 that not only mimics the corresponding row of G but also accounts for the value of the fourth entry of X, here 1. In general, with perfect phase and/or amplitude control, this additional encoding of X is trivial. Yet, realistic scenarios often imply rather limited control, such as phase-only. Then, since the reflection coefficient of the wavefront shaper elements cannot trivially be set to zero, a configuration creating destructive interferences has to be found; this is why the schematic does not simply display the segments corresponding to zero entries of X as zero.
To demonstrate the wide applicability to all wave types and the ease of implementing our proposed concept, we consider the original case of using everyday wireless communication signals, such as Wi-Fi waves, bouncing around inside an indoor room. Indeed, indoor environments constitute cavities of low quality factor for Wi-Fi waves. Their geometries are in general not perfectly regular, such that these cavities may be labelled "disordered" from the wave's perspective. The reflections off the irregular cavity walls interfere and give rise to a speckle-like wave field [29], very much like scattering events in multiply scattering optical media such as thin paint layers or biological tissue.
A convenient way of probing the field in a cavity is to measure the Green's function between two independent positions [30, 31]. We need as many independent single-frequency transmission measurements as Y is supposed to have entries. In principle, placing antennas at positions that are at least half a wavelength apart is an obvious option; one may also use a single pair of antennas and working frequencies that are sufficiently different [32]. In the realistic scenario of Wi-Fi, the channel state information (CSI) between an access point and a wireless device is exactly what we need: a complex-valued transmission measurement between two points inside the cavity. Modern multi-user telecommunication networks routinely sound channels, for instance using beacon signals, for downlink beamforming [33]. CSI can be accessed with off-the-shelf Intel 5300 cards, and has successfully been extracted and used for instance as fingerprint for object localization [34].
The crucial remaining question is: How can a Wi-Fi wave field in an indoor environment be shaped? The concept of wavefront shaping originally emerged in optics but has recently been transposed to the microwave domain using simple phase-binary reflect-arrays, with important applications in telecommunication, imaging and energy transfer [35][36][37]. Each element of such an array can be configured electronically to mimic a Dirichlet or Neumann boundary condition. Our reflectarray is a home-made tunable metasurface whose working principle is based on the hybridization of resonances (see Ref. [38] and Methods for further details). Stated differently, we can program the phase shift of the reflected wave for each of the 96 elements individually to be 0 or π [39,40].
In our laboratory implementation, in order to work under stable, well-controlled conditions, we mimic the indoor room using a disordered metallic cavity with numerous absorbers placed on its walls such that its quality factor Q = 179 is comparable to the one we measured in an office room (see Methods). A schematic is shown in Fig. 2. Moreover, we mimic CSI between wireless devices and a base station by transmission measurements between independent monopole antennas with a vector network analyzer.
We seek to perform the following 4 × 4 operation: (1) While we could use any complex-valued matrix for our demonstration, we deliberately selected the one in Eq. 1 because it describes the (discrete) Fourier transform operation. Incidentally, complex media in conjunction with WFS in optics, or Time Reversal in acoustics, have repeatedly been described as "opaque lens" or "scattering lens" in the literature [14,41], alluding to the most basic property of a traditional lens: it focuses light in geometrical optics. In our present work, we use a complex medium to Fourier transform an impinging wavefront, such that it becomes possible to extend the analogy from the ray picture to the properties of a traditional lens in the realm of "wave" optics. While a lens in optics is a standard component, transposed to the microwave domain our concept unveils yet more of its intriguing properties: (i) lenses are not as widely-used; (ii) reflectarray and antennas, mimicking input and output of the Fourier transform, are not aligned in any way but placed randomly and the medium carrying out the operation is around rather than between inputs and outputs. Finally, we note that Eq. 1 also corresponds to a 2-qubit Quantum-Fourier-Transform (QFT), the pivotal step in Shor's famous factoring algorithm [42] -although of course our experiment is entirely classical. Much effort to experimentally calculate the results of a QFT has been reported with cumbersome setups involving, for instance, trapped ions or nuclear magnetic resonance, to name a few [43,44].
In a one-off initiation step, we first characterise the impact of each reflect-array element on each measured transmission pair. This yields a matrix H that we term Impact Matrix (IM) to avoid ambiguities. In transmission geometries, such as in our schematic in Fig. 1 and in standard optical WFS experiments, the IM is simply the medium's Transmission Matrix [17]. The Physics of our experimental implementation in a microwave cavity (Fig. 2) is however richer than captured by a simple, linear transmission formalism. Since the reverberating wave revisits the reflect-array multiple times, the impact of a given element on the wave field is correlated to some extent to the configuration of the remaining reflect-array elements: H is itself a function of the reflect-array configuration. A reflect-array element acts not like a single but multiple (secondary) sources, due to reverberation.
Even in cavities of low quality factor, such as an indoor environment, these long-range correlations remain substantial (see SI Fig. S1). For the sake of mathematical and experimental simplicity, we measure a first-order approximation of H that ignores long-range correlations (see Methods). Then, aware of the limited predictive power of such an approximative IM, we take appropriate measures to overcome the resulting imperfections. First, we numerically identify a library of segment configurations that appropriately tame the wavefront as aforementioned (see Methods for the detailed procedure), supposing our first-order IM is exact. Of course, the computation outcomes that we expect based on our IM, Y pred , are not exactly what we will measure, Y meas , since we did not account for the long-range correlations yet. The effect of the latter can be interpreted as a random phasor For a given realization, δ is a deterministic but seemingly arbitrary random complex number that results in an inevitable inaccuracy of our computation outcome Y meas . A realization may be defined as one instance of the experiment with some specific (disordered) cavity geometry. A different realization will have the same global parameters (cavity volume, quality factor, . . . ) but a different geometry and thus a different random δ. Therefore, we can exploit the realization-dependence of δ to average out the intrinsic inaccuracy: δ realizations = 0, and thus Y meas realizations = Y pred . An illustration of this effect based on experimental data is provided in the Supplementary Information. By averaging over realizations of disorder we are thus able to mitigate the adverse effects of using a first-order IM.
Realizations of disorder are easily achieved in our experiment by rotating a mode-stirrer by 12 • (see Fig. 2) or changing the operating frequency by more than a correlation frequency (see Methods). Both concepts could of course be used in an analogous manner with Wi-Fi in an indoor room. The constant motion of inhabitants naturally realizes disorder. A single realization would have to be carried out within the coherence time of the medium, which is envisageable in sight of MHz rates at which the PIN-diodes controlling the reflect-array elements can be switched with improved electronics [45]. Alternatively, using the CSI between different pairs of access point and mobile device (e.g. if there are many mobile devices in a room) is another way to realize disorder without relying on motion inside the room.
En passant, we note that averaging over realizations of disorder is also extremely useful in case of very limited wavefront control. For instance, each segment in our 4 × 4 experiment is made up of only 24 phase-binary elements which is insufficient to perfectly tame the random walk. For a given realization, even the predicted computation outcomes, Y pred , based on the IM such that the long-range correlation issue does not intervene, do not perfectly match the theoretically expected computation outcomes, Y exact , due to the limited control over the wavefront. This difference, once again, can be interpreted as a random realization-dependent phasor that is averaged out by realizing disorder. Hence, the idea of averaging over realizations of disorder serves two purposes in our case.
In Fig. 3(a-f) we present our experimentally obtained computation results Y for six different input vectors X, based on 150 realizations of disorder (30 mode-stirrer positions, 5 independent working frequencies). We observe an excellent agreement with the theoretically expected results, the average computation error = |Y exact j −Y j | j being 3.8%, averaged over all possible input vectors X. Essentially, by configuring the reflect-array according to the previously established library for a given input X, we create a cavity geometry in which the transmission between the selected antenna pairs corresponds to the computation output Y . The dependence of on the number of realizations displayed in Fig. 3(g) confirms that the inaccuracy due to the long-range correlations can indeed be averaged out successfully by realizing disorder. The same is true for the issue related to the very limited control over the wavefront, as seen with the purple curve solely based on the IM: is reduced from 8.3% to 0.7%.
Having demonstrated experimentally that Wi-Fi waves in an indoor environment could perform a 4 × 4 computational operation, we now address the question whether, given our system, we can perform even larger operations. If we were to bluntly apply our previous procedure, identifying appropriate configurations for each segment would become doubly more difficult: (i) given the fixed reflect-array size, each wavefront segment would consist of even fewer elements, further reducing the available degrees of freedom; (ii) each segment configuration would have to satisfy even more constraints, since there are also more observation points due to the larger size of the output vector Y . Averaging over more realizations of disorder could cushion these adverse effects to some extent, but there is a more elegant solution. We can exploit the linearity of the computation operation Y = GX in which the entries of Y are independent from each other. Thereby, we compute one entry of Y after the other, such that we convert one M × N operation into M 1 × N operations (see SI for illustration) that are performed one after the other. This is possible thanks to the reconfigurability of our wave-based analogue computation unit -a fundamental advantage of our proposed system over metamaterial computation units fabricated for a single specific task. We present in Fig. 4 the experimentally obtained Fourier transform of the silhouette of two Parisian monuments. We We have demonstrated experimentally that any disordered medium can be employed as reconfigurable analogue computation unit, subject to appropriate wavefront shaping -even ubiquitous Wi-Fi waves can perform analogue computation simply by reverberating inside an indoor room.

ADDITIONAL INFORMATION Data Availability
The raw data is available at ZENODO-DOI. Codes may be obtained from the corresponding author upon reasonable request.

Competing Interests
The authors declare no competing interests.

METHODS
Experimental Setup. We use a metallic cavity (V = 1.01 × 0.86 × 1.28 m 3 = 1.1 m 3 ) that is disordered due to a few irregularly shaped metallic objects glued to the walls as well as the presence of the large mode-stirrer. Pieces of electromagnetic absorbers (10 × 10 × 5 cm 3 ) are isotropically distributed across the cavity's walls to lower its quality factor from 1023 to 179. This is on the order of the quality factor of 234 that we have measured in an office room in our laboratory. We estimate Q as πf 0 /µ, where µ is the exponential decay constant extracted from the Green's function envelope averaged over 30 mode-stirrer positions. Our cavity's correlation frequency is thus ∆f corr = f 0 /Q ≈ 15 MHz. Consequently we choose five working frequencies (within the reflect-array's operating band) that are separated by at least ∆f corr in order for them to constitute independent realizations. The reflect-array consists of n = 96 elements whose design is essentially the one presented in Ref. [38].
However, each element has independent control over the two polarizations of the electromagnetic field, rather than acting only on one polarization as in Ref. [38]. We synchronize each element's horizontal and vertical polarization configuration to simplify our experiment. Each of our monopole antennas picks up the field polarized along its axis; since polarizations are mixed in a chaotic cavity we have n = 96 controllable elements whose average impact on the wave field is twice as much as it would be if they controlled only one polarization.
Measuring the Impact Matrix. The entry H i,j quantifies the impact of the reflect-array element indexed i on the transmission between the pair of antennas indexed j. We measure the Impact Matrix (IM) using the Hadamard rather than the canonical basis, following the procedure used in Ref. [17]. The elements of the Hadamard basis being +1 or −1, this is a perfect match to our experiment where we can access exactly these two values for each element. Moreover, working in the Hadamard basis has a favourable averaging effect as illustrated in the SI Fig. S1 and benefits furthermore from an improved signal-to-noise ratio. One important difference to Ref. [17], however, is that we only control a portion of the wave field with our reflect-array covering about 6% of the cavity surface. The cavity wave field at the working frequency f 0 may be decomposed as sum of the cavity's N = 1 Q ≈ 96 modes that overlap at the working frequency f 0 due to their finite line-widths, c being the speed of light. As discussed in Ref. [47], the reflect-array's control over the wave field may be modeled as equivalent to effectively controlling some p = 2 × (n/3) ≈ 64 modes, where n = 96 is the number of reflect-array elements and the factor of 2 accounts for the independent control over both polarizations. Thus, given that we effectively only control about p/N ≈ 67% of the cavity modes, to successfully define an IM we subtract from each transmission Y the static cavity contribution U that cannot be altered. We obtain U as Y random V , where V is a vector defining the reflect-array configuration, and then we work with Y = HX, where Y = Y − U as in Ref. [28]. Note that this subtraction is in principle avoidable by using a larger reflect-array such that p ≈ N ; then, the entire wave field is modulated.
Optimisation Algorithm. Given the IM, we numerically identify the configurations that yield the desired tamed walks illustrated in Fig. 1(b). Consider the 4 × 4 example given in Fig. 1. Let V D be the part of vector V that belongs to segment D, let Y i D denote the contribution of group D to observation point i, and let h i D be the corresponding part of the IM H that links V D to Y i D . Note that selecting the members of segment D amongst all wavefront shaper elements does not have to be done in a regular manner as in our schematic in Fig. 1. In our experiment, we make the selection randomly; a careful selection could in fact be a further parameter of the optimisation. The objective is to find a V D , given h 1 D to h 4 D , such that simultaneously specific values of Y 1 D to Y 4 D are obtained. This can be expressed as a single stacked equation, In the case of perfect phase and amplitude control over each wavefront element, the equation is easily In our experimental setup, however, the entries of V may only take the values of +1 or −1 due to our phase-binary only control over the wavefront. A plethora of methods to solve this ill-posed problem exists [48,49]. Here we opt for a home-made solver inspired by standard iterative sequential optimization algorithms used in optical wavefront shaping [50]. First, we define a normalization factor γ between the entries of G and the target amplitudes of the tamed walks. We choose γ = m 1/2 |H i,j | i,j , where m is the number of elements in a given wavefront segment. This choice corresponds to the expected value of the resultant of a random walk of m steps of step size |H i,j | i,j . Next, we define a cost function CF as distance in the complex plane between an entry of γG D,j and the corresponding value of h j D V D , averaged over all output points j. Starting with a random configuration of V D , we test element by element if flipping its state reduces CF in which case we update V D accordingly. We run 20 loops over all elements. Finally, we repeat the procedure with 250 different random initial V D and select the overall best final V D yielding the lowest CF as entry for the library.
As outlined in the main text, we opt for the scenario of additionally encoding the impinging wavefront into the wavefront shaper configuration. Since we do not have perfect phase and amplitude control on the reflect-array, this is not trivial. To encode an input vector entry α other than unity into, for example V D , we identify a configuration for V D such that its contribution h j D V D to the observation point indexed j is αγG D,j , for all j. For Fig. 3, we restricted ourselves to α ∈ {0, 1}; for  About 6% of a disordered metallic cavity's surface are covered with a reflect-array consisting of 96 pixels with tunable reflection coefficient (r = ±1); see Ref. [38] and Methods for detailed working principle.
Numerous electromagnetic absorbers placed on the walls lower the cavity quality factor to Q = 177, comparable to an indoor room. With a vector network analyzer, we measure the transmission between four simple monopole Wi-Fi antennas, mimicking channel state information between an access point and four mobile devices. A mode-stirrer, rotated in steps of 12 • , enables realizations of disorder.   For the interested reader, here we provide numerous additional details that complement the manuscript and provide further illustrations. This document is organized as follows: A. Long-range correlations in the Impact Matrix.
B. Averaging over Realizations of Disorder.

C. Sequential Analogue Computation.
A To now obtain many measurements of H i,j , we work in the usual, canoncial basis. We define V + and V − as being equal to the same random phase-binary configuration V , except for the ith entry that is fixed to be +1 or −1 for V + and V − , respectively. By measuring the corresponding transmissions Y + and Y − , we can obtain an estimate of H i,j as follows: H can i,j = (Y + − Y − )/2. We repeat this 160 times with different random vectors V , yielding 160 estimates of H can i,j . From those, we can further define an average over the 160 estimates: H av i,j = H can i,j . In Fig. S1 we display for two examples in the complex plane the experimentally measured cloud of 160 points for H can i,j , their average H av i,j as well as H had i,j . Note that the reflect-array element considered in (b) is defect, thus certainly having no impact on the cavity wave field.
The long-range correlations are clearly visualized by the non-negligible radius of the blue cloud.
The numerical identification of the appropriate configuration for a given wavefront segment works with one value of H i,j . Yet, the actual value in the computation phase will depend on how the reflect-array is configured then. Clearly, using H av i,j in the numerical optimisation is the best choice. It is, however, very cumbersome to measure. The measurement in the Hadamard basis that we opt for requires only a single measurement per IM entry and has an averaging effect because it employs all reflect-array elements. It therefore seems to be a good trade-off between measurement duration and accuracy. The averaging over realizations is inevitable in any case, but the further from H av i,j the value used in the numerical optimisation is, the more likely it is that more realizations over disorder are needed to converge to the exact computation results.
Finally, we note that in Ref. [1] a first-order approximation of an IM in a similar system (low-Q disordered microwave cavity) was measured and successfully used for focusing without averaging over realizations of disorder. This worked because said focusing experiment only required a phasebinary decision in order to align each element with the overall desired direction in the complex plane. As seen in Fig. S1(a), the cloud usually remains within one quadrant, here the upper right one. The focusing decision thus only requires a choice as to whether we want some value in the upper right (element configuration r = +1) or lower left (element configuration r = −1) quadrant to maximize the length of the resultant. In the present work, the optimisation constraints are significantly more challenging to satisfy: a specific amplitude and phase has to be achieved with the tamed walk, simultaneously at several observation points.

B. Averaging over Realizations of Disorder
In Fig. S2, for the experiment of computing the Fourier transform of the Eiffel Tower's silhouette (see main text Fig. 4(a)), we display for each output point the outcomes of the 150 individual realizations (30 mode-stirrer positions, 5 independent working frequencies), as well as their average and the exact value that are shown in Fig. 4(a).  Fig. 4(a) in the main text.

C. Sequential Analogue Computation
In this section, we first visually illustrate the principle of sequential analogue wave-based computation explained in the main text. Then, we provide the 16 × 16 complex-valued operation G that we implemented, corresponding to Fig. 4 in the main text. Finally, we show some further computation examples obtained with said operation that we did not include in the main text for conciseness.

Visual illustration of Sequential Analogue Computation
Here, we provide a visualization of how we can achieve a much bigger effective computation unit than our fixed-size physical system. Consider the general matrix formalism for an operation G as shown in Fig. S3(a). For an M × N operation G, each output entry in Y is in fact the product of the corresponding row of G, itself a vector, and the input vector X. We can thus separately compute each output entry in a 1 × N operation. Since the resultants of the complex walks (main text Fig. 1) have to be tamed only for one observation point at any given time, the optimization constraints are considerably relaxed such that it becomes possible to implement larger operations with more inputs and outputs than in a single-shot computation scheme. The procedure is thus essentially looped over each of the output vector entries, as schematically shown in Fig. S3(b). Moreover, this means that we only need one transmission measurement for the sequential version. Note that one may in principle also compute more than one of the output vector's entries per loop to speed up the computation time.
terms of the corresponding intensities (right). The results were obtained over 150 realizations of disorder like the results in the main text.