Methods for Increasing Transmission Through Scattering Random Media

Abstract

The method, system, and computer-readable medium increase transmission of waves through a highly scattering random medium. Transmission is increased by iteratively refining wavefronts using measurements of the backscatter wavefronts resulting from transmission of waves into the medium. The process of double phase conjugation by time-reversing a wavefront, transmitting the time-reversed wavefront into the medium, and time-reversing the reversing backscatter wavefront is leveraged to implement the method in a physical system using a phase conjugate mirror. In an embodiment, transmission may be increased by phase-only modulation of the wavefronts. In an embodiment, the invention may be used to focus transmission through the medium to a location opposite the wave source.

Description

TECHNICAL FIELD

The invention relates generally to methods and devices for signal transmission and, more particularly, to methods and devices for transmitting waves through highly scattering random media.

BACKGROUND

Transparent media such as glass and air are transparent because light propagates through them without being scattered or absorbed. In contrast, opaque materials prevent transmission of light by absorption or scattering. Highly scattering random media such as biological tissue, turbid water, white paint, and egg shells are opaque because they contain randomly arranged particles, which cause light to scatter in many directions. This scattering effect becomes more pronounced as the thickness of the highly scattering random medium increases, and less and less normally incident light is transmitted through the medium. Scattering thus severely curtails the usefulness of optical techniques for sensing and imaging in highly scattering random media.

Previous work in the field has shown the existence of eigen-wavefronts that dramatically increase the transmission of electromagnetic signals in highly scattering random media (viz., wavefronts with transmission coefficients close to 1). By measuring the scattering of light that has passed through a highly scattering random medium, such eigen-wavefronts can be determined. For a given highly scattering random medium, an optimal wavefront that maximizes transmission through the medium can be produced from the eigen-wavefronts. Existing coordinate descent methods for determining such optimal wavefronts suffer from two main problems that limit the use of such methods in many applications. First, existing methods require measurement within or on the far side of the medium from the source of the signal. Second, existing methods converge to an optimal wavefront slowly. Existing coordinate descent methods maximize the measured intensity of each mode of a transmitted wavefront while holding the amplitudes and phases of all other modes constant. The process must then be repeated for each of M modes, requiring on the order of M measurements. Other existing methods use multiple frequencies to find the optimal phases simultaneously, but these other methods still require measurement through the medium and have only been shown to work for small numbers of modes.

SUMMARY

Disclosed is a method, system, and computer-readable medium storing instructions for increasing transmission of a wavefront through a highly scattering random medium using measurements of the backscatter field produced by wavefronts incident upon the medium. One embodiment consists of a computer-implemented method including refining a test wavefront in one or more iterations to produce a high-transmission wavefront, terminating the iterative refinement upon reaching a numerically prespecified level of transmission, determining a modified wavefront to transmit through the medium, and transmitting the modified wavefront through the medium. Another embodiment consists of a system including a processor, a spatial light modulator configured to modulate the wavefront and transmit the modulated wavefront into the medium, a detector configured to measure the backscatter field resulting from the transmission of the wavefront into the medium, and a program memory storing executable instructions that when executed by the one or more processors cause the system to refine the wavefront in one or more iterations, terminate the iterative refinement upon reaching a numerically prespecified level of transmission, modify the wavefront based upon the final test wavefront, and transmit the modified wavefront through the medium. Another embodiment consists of a tangible, non-transitory computer-readable medium storing instructions that when executed by one or more processors of a computer system cause the computer system to refine the wavefront in one or more iterations, terminate the iterative refinement upon reaching a numerically prespecified level of transmission, modify the wavefront based upon the final test wavefront, and transmit the modified wavefront through the medium.

The iterative refinement of the wavefront includes transmitting a wavefront into the medium, receiving a measurement of the back-scatter, determining whether the level of transmission through the medium has reached a numerically prespecified level to terminate the iterative refinement, and refining the wavefront for transmission into the medium in the next iteration. In an embodiment, the iterative refinement of the wavefront may consist of transmitting a first test wavefront into the medium, receiving a measurement of the first backscatter wavefront resulting from the first test wavefront, time-reversing the first backscatter wavefront to obtain a second test wavefront, transmitting the second test wavefront into the medium, receiving a measurement of the second backscatter wavefront resulting from the first test wavefront, time-reversing the second backscatter wavefront to obtain a third test wavefront, determining whether the level of transmission through the medium of the third test wavefront has reached a numerically prespecified level to terminate the iterative refinement, and refining the wavefront for transmission into the medium in the next iteration using the third test wavefront. In some embodiments, a phase conjugate mirror is used to time-reversing the first and second backscatter wavefronts. In an embodiment, the second test wavefront may be determined using only the phase of the modal coefficients of the first backscatter wavefront. Some embodiments may use phase and amplitude modulation of the wavefronts or may use phase-only modulation.

Refining the test wavefront in each iteration may involve adjusting the test wavefront in proportion to the third test wavefront. Refining the test wavefront may also involve adjusting the test wavefront in proportion to the average magnitude of the modal coefficients of the first backscatter wavefront. Refining the test wavefront may also involved adjusting the phase of the test wavefront in proportion to the average magnitude of the modal coefficients of the first backscatter wavefront. In an embodiment, refining the test wavefront may further comprise determining a residual vector using the third test wavefront.

In an embodiment, the refinement may include determining a matrix H using the first and third test wavefronts and a matrix Q, the columns of which are the modal coefficient vectors of the test wavefronts of each of the one or more iterations. The eigenvectors of the matrix H and the matrix Q may be used to modify the wavefront to focus the wavefront through the medium.

BRIEF DESCRIPTION OF THE DRAWINGS

The figures described below depict various aspects of the applications, methods, and systems disclosed herein. It should be understood that each figure depicts an embodiment of a particular aspect of the disclosed applications, systems and methods, and that each of the figures is intended to accord with a possible embodiment thereof. Furthermore, wherever possible, the following description refers to the reference numerals included in the following figures, in which features depicted in multiple figures are designated with consistent reference numerals.

FIG. 1 illustrates an example highly scattering random medium.

FIG. 2 illustrates an example system for increasing transmission through a highly scattering random medium.

FIG. 3 illustrates an example computer system for implementing the methods for increasing transmission through a highly scattering random medium.

DETAILED DESCRIPTION

Although the following text sets forth a detailed description of numerous different embodiments, it should be understood that the legal scope of the invention is defined by the words of the claims set forth at the end of this patent. The detailed description is to be construed as exemplary only and does not describe every possible embodiment, as describing every possible embodiment would be impractical, if not impossible. One could implement numerous alternate embodiments, using either current technology or technology developed after the filing date of this patent, which would still fall within the scope of the claims.

It should also be understood that, unless a term is expressly defined in this patent using the sentence “As used herein, the term ‘______’ is hereby defined to mean . . . ” or a similar sentence, there is no intent to limit the meaning of that term, either expressly or by implication, beyond its plain or ordinary meaning, and such term should not be interpreted to be limited in scope based on any statement made in any section of this patent (other than the language of the claims). To the extent that any term recited in the claims at the end of this patent is referred to in this patent in a manner consistent with a single meaning, that is done for sake of clarity only so as to not confuse the reader, and it is not intended that such claim term be limited, by implication or otherwise, to that single meaning. Finally, unless a claim element is defined by reciting the word “means” and a function without the recital of any structure, it is not intended that the scope of any claim element be interpreted based on the application of 35 U.S.C. §112(f).

Although the description within this application discusses transmission in terms of light or optics, it should be understood that these terms are used for convenience and brevity. The method and system apply equally to any electromagnetic, acoustic, optical, or other waves propagating through a highly scattering random medium. As used herein, the term “highly scattering random medium” means any medium that has a low transmission coefficient to most waves in the applicable range of frequencies, that causes a high degree of scattering within the medium in the applicable range of frequencies, and that has a low degree of absorption within the medium in the applicable range of frequencies. The applicable range of frequencies may depend upon the particular application to which the described method or system is applied. As used herein, the term “medium” refers to a transmission medium. As used herein the superscripts ^H and * applied to a matrix or vector (e.g., X^H or X*) denote the complex conjugate transpose and the complex conjugate, respectively, of the matrix or vector to which they refer.

Transmission in Highly Scattering Random Media

FIG. 1 illustrates an example highly scattering random medium 100 as a two-dimensional periodic slab of thickness L and periodicity D. The slab's unit cell occupies the space 0≦x<D and 0≦y<L and contains N, infinite and z-invariant circular cylinders Cyl that are placed randomly within the cell. The circular cylinders may be either perfect electrically conducting or dielectric with refractive index n_d. In the example medium, the cylinders do not overlap and all fields are TM_z polarized. Electric fields in the halfspace where y<0 (i=1) and the halfspace where y>L (i=2) are denoted as e_i(ρ)=e_i (ρ){circumflex over (z)}. These fields' complex amplitudes e_i(ρ) can be decomposed in terms of propagating in the +y and −y directions as e_i(ρ)=e_i⁺(ρ)+e_i⁻(ρ), where

$\begin{array}{cc}{e}_i^{\pm }\left(\underset{\_}{\rho }\right)=\sum _{n=-N}^Nh_na_{i,n}^{\pm }{}^{-j{\underset{\_}{k}}_n^{\pm }·\underset{\_}{\rho }}.& \mathrm{Equation}1\end{array}$

The components of Equation 1 above are defined as follows: ρ=x{circumflex over (x)}+yŷ≡(x,y), k_n^±=k_n,x{circumflex over (x)}±k_n,yŷ≡(k_n,x,±k_n,y), k_n,x=2πn/D, k_n,y=2π√{square root over ((1/λ)²−(n/D)²)}{square root over ((1/λ)²−(n/D)²)}, λ is the wavelength, and h_n=√{square root over (∥k_n^±∥₂/k_n,y)} is a power-normalizing coefficient. It is assumed that N=[D/λ], so only propagating waves are modeled. The number of modes of propagating waves is denoted as M=2N+1, and the modal coefficients are denoted a_i,n^± where i=1, 2 and n=N, . . . , 0, . . . , N. Vectors of the modal coefficients are denoted a_i^±=[a_i,−N^± . . . a_i,0^± . . . a_i,N^±]^T where T denotes transposition. The modal coefficients a_i,n⁺ are related by the scattering matrix S as follows:

$\begin{array}{cc}\left[\begin{array}{c}{\underset{\_}{a}}_1^{-}\\ {\underset{\_}{a}}_2^{+}\end{array}\right]=\underset{\underset{=:S}{}}{\left[\begin{array}{cc}{S}_{11}& S_{12}\\ S_{21}& S_{22}\end{array}\right]}\left[\begin{array}{c}{\underset{\_}{a}}_1^{+}\\ {\underset{\_}{a}}_2^{-}\end{array}\right],& \mathrm{Equation}2\end{array}$

For the present disclosure, the highly scattering random medium 100 is excited only from the direction of y<0; hence, a₂⁻=0 and S₁₂ and S₂₂ may be disregarded. Thus, the modal coefficients of the backscatter wavefront are a₁⁻=S₁₁·a₁⁺, and the modal coefficients of the transmitted wavefront are a₂⁺=S₂₁·a₁⁺. For any incident field amplitude e₁⁺(ρ), the transmission coefficient may be defined as

$\begin{array}{cc}\tau \left({\underset{\_}{a}}_1^{+}\right):=\frac{{S_{21}·{\underset{\_}{a}}_1^{+}}_2^2}{{{\underset{\_}{a}}_1^{+}}_2^2},& \mathrm{Equation}3\end{array}$

and the reflection coefficient may be defined as

$\begin{array}{cc}\Gamma \left({\underset{\_}{a}}_1^{+}\right):=\frac{{S_{11}·{\underset{\_}{a}}_1^{+}}_2^2}{{{\underset{\_}{a}}_1^{+}}_2^2},& \mathrm{Equation}4\end{array}$

The transmission coefficient of a normally incident wavefront may be denoted by T_normal=τ([0 . . . 0 1 0 . . . 0]^T).

From the foregoing, it can be seen that the maximum transmission of power through a highly scattering random medium 100 may be obtained by designing an optimal wavefront a_opt incident to the highly scattering random medium 100. The optimal wavefront a_opt is a wavefront for which the transmission coefficient obtains its maximum value in a system with scattering matrix S. As such, the optimal wavefront may be determined as follows:

$\begin{array}{cc}\begin{array}{c}{\underset{\_}{a}}_{\mathrm{opt}}=\underset{{\underset{\_}{a}}_1^{+}}{\mathrm{argmax}}\tau \left({\underset{\_}{a}}_1^{+}\right)\\ =\underset{{\underset{\_}{a}}_1^{+}}{\mathrm{argmax}}\frac{{S_{21}·{\underset{\_}{a}}_1^{+}}_2^2}{{{\underset{\_}{a}}_1^{+}}_2^2}\\ =\underset{{{\underset{\_}{a}}_1^{+}}_2=1}{\mathrm{argmax}}{S_{21}·{\underset{\_}{a}}_1^{+}}_2^2\end{array}& \mathrm{Equation}5\end{array}$

where ∥a₁⁺∥₂=1 represents the incident power constraint.

The highly scattering random medium 100 can be approximated by a lossless medium. In a lossless medium, the scattering matrix S will be unitary; i.e. S^H·S=I where I is the identity matrix. Consequently, S₁₁^H·S₁₁+S₂₁^H·S₂₁=I, and the optimal wavefront may also be determined as follows:

$\begin{array}{cc}\begin{array}{c}{\underset{\_}{a}}_{\mathrm{opt}}=\underset{{{\underset{\_}{a}}_1^{+}}_2=1}{\mathrm{argmax}}\underset{\underset{={\left({\underset{\_}{a}}_1^{+}\right)}^H·\left(I·S_{11}^H·S_{11}\right)·{\underset{\_}{a}}_1^{+}}{}}{{\left({\underset{\_}{a}}_1^{+}\right)}^H·S_{21}^H·S_{21}·{\underset{\_}{a}}_1^{+}}\\ =\underset{{{\underset{\_}{a}}_1^{+}}_2=1}{\mathrm{argmax}}{S_{11}·{\underset{\_}{a}}_1^{+}}_2^2\\ =\underset{{\underset{\_}{a}}_1^{+}}{\mathrm{argmax}}\Gamma \left({\underset{\_}{a}}_1^{+}\right).\end{array}& \mathrm{Equation}6\end{array}$

In a lossless medium, the wavefront that minimizes backscattering also maximizes transmission. The same is true in a highly random scattering medium because it approximates a lossless medium. The present invention uses this equivalence to develop a method and system of maximizing transmission by minimizing backscatter, which enables the determination of the optimal wavefront a_opt using only the backscatter wavefront with modal coefficient vector a₁⁻.

Additionally, or alternatively, the optimal wavefront a_opt may be determined using phase modulation without amplitude modulation. Such an approach is of particular benefit in applications where amplitude modulation may be impractical, such as applications in the optical range, but the technique may be applied more generally. The modal coefficient vector a₁⁺ can be separated into amplitude and phase components, where |a_1,n⁺| and ∠a_1,n⁺ denote the magnitude and phase of modal coefficient a_1,n⁺=|a_1,n⁺| exp(j∠a_1,n⁺), respectively. The phase-only modal coefficient vector may then be defined as ∠a₁⁺=“∠a_1,−N⁺ . . . ∠a_1,0⁺ . . . ∠a_1,N⁺]^T. A set P_c^M may be defined comprising phase vectors each with M elements and non-zero positive real-valued constant c of the form

$\begin{array}{cc}\underset{\_}{p}\left(\underset{\_}{\theta };c\right)={\sqrt{\frac{c}{M}}\left[{}^{{\mathrm{j\theta }}_{-N}}\dots {}^{{\mathrm{j\theta }}_0}\dots {}^{{\mathrm{j\theta }}_N}\right]}^T,& \mathrm{Equation}7\end{array}$

where θ=[θ_−N . . . θ₀ . . . θ_N]^T is an M-vector of phases. With the incident power constraint ∥a₁⁺∥₂=1 as above, the optimal wavefront a_opt can be determined using phase-only modulation as

$\begin{array}{cc}{\underset{\_}{a}}_{\mathrm{opt}}=\underset{{\underset{\_}{a}}_1^{+}\in P_c^M}{\mathrm{argmax}}\tau \left({\underset{\_}{a}}_1^{+}\right)=\underset{{\underset{\_}{a}}_1^{+}\in P_c^M}{\mathrm{argmax}}\frac{{S_{21}·{\underset{\_}{a}}_1^{+}}_2^2}{{{\underset{\_}{a}}_1^{+}}_2^2}=\underset{{\underset{\_}{a}}_1^{+}\in P_1^M}{\mathrm{argmax}}{S_{21}·{\underset{\_}{a}}_1^{+}}_2^2.& \mathrm{Equation}8\end{array}$

The incident power constraint requires that c=1 and, therefore, a₁⁺∈P₁^M and a_opt∈P₁^M. Where c=1, the additional phase vector form p(θ):=p(θ; 1) may also be defined. The optimal wavefront may then be determined as a phase-only modulated wavefront a_opt=p(θ_opt), where the optimal phase vector θ_opt may be determined as a phase vector θ for which transmission coefficient achieves its maximum value for a highly scattering random medium 100 with scattering matrix S, as

$\begin{array}{cc}{\underset{\_}{\theta }}_{\mathrm{opt}}=\underset{\underset{\_}{\theta }}{\mathrm{argmax}}\underset{\underset{={\left(\underset{\_}{p}\left(\underset{\_}{\theta }\right)\right)}^H·\left(1-S_{11}^H·S_{11}\right)\underset{\_}{p}\left(\underset{\_}{\theta }\right)}{}}{{\left(\underset{\_}{p}\left(\underset{\_}{\theta }\right)\right)}^H·S_{21}^H·S_{21}·\underset{\_}{p}\left(\underset{\_}{\theta }\right)}=\underset{\underset{\_}{\theta }}{\mathrm{argmin}}{S_{11}·\underset{\_}{p}\left(\underset{\_}{\theta }\right)}_2^2=\underset{\underset{\_}{\theta }}{\mathrm{argmin}}\Gamma \left(\underset{\_}{p}\left(\underset{\_}{\theta }\right)\right).& \mathrm{Equation}9\end{array}$

As above, the highly scattering random medium 100 can be approximated as a lossless medium, in which case the optimal phase vector θ_opt may be obtained instead by minimizing the backscatter (represented by the reflection coefficient Γ) as

$\begin{array}{cc}{\underset{\_}{a}}_{\mathrm{opt}}=\underset{{\underset{\_}{a}}_1^{+}\in P_1^M}{\mathrm{argmin}}{S_{11}·{\underset{\_}{a}}_1^{+}}_2^2,& \mathrm{Equation}10\end{array}$

In a lossless medium, using amplitude and phase modulation as discussed above is expected to yield an average optimal transmission efficiency as measured by the transmission coefficient τ of approximately 1. Due to the limitation on amplitude, phase-only modulation may not achieve perfect transmission, but phase-only modulation as described above is expected to yield an average optimal transmission efficiency as measured by the transmission coefficient τ of at least approximately π/4 in a lossless medium.

Iterative Algorithm for Maximizing Transmission Using Backscatter Analysis

The optimal wavefront a_opt or the optimal phase vector θ_opt can be determined in a physically realizable iterative system using backscatter analysis, thereby eliminating the necessity of measuring the intensity of the transmitted wavefront that has passed through the highly scattering random medium 100 to the halfspace where y>L. In each iteration of the several algorithms below, the estimation of the optimal wavefront a_opt or the optimal phase vector θ_opt is directly or indirectly updated until a numerically prespecified level of transmission is reached. Several embodiments of the system and method are discussed below.

One method of determining the optimal wavefront a_opt involves a steepest descent algorithm as set forth in Table 1 below. The steepest descent algorithm utilizes the negative gradient of the objective function ∥S₁₁·a₁⁺∥₂² from Equation 10 as a search direction to converge upon the optimal wavefront a_opt by iterative steps of predetermined step size μ>0. For any kth iteration of the steepest descent algorithm, the next iteration k+1 is determined by a step of size μ in the direction of the negative gradient of the objective function, expressed as

$\begin{array}{cc}\begin{array}{c}{\underset{\_}{\tilde{a}}}_{1,\left(k\right)}^{+}={\underset{\_}{a}}_{1,\left(k\right)}^{+}-\mu \frac{\partial {S_{11}·{\underset{\_}{a}}_1^{+}}_2^2}{\partial {\underset{\_}{a}}_1^{+}}{|}_{{\underset{\_}{a}}_1^{+}={\underset{\_}{a}}_{1,\left(k\right)}^{+}}\\ ={\underset{\_}{a}}_{1,\left(k\right)}^{+}-2\mu S_{11}^H·S_{11}·{\underset{\_}{a}}_{1,\left(k\right)}^{+},\end{array}& \mathrm{Equation}11\end{array}$

where a_1,(k)⁺ represents the modal coefficient vector of the wavefront produced at the kth iteration of the algorithm. In each iteration, a refined wavefront with modal coefficient vector ã_1,(k)⁺ is determined as indicated in line 6 of Table 1. The refined wavefront with modal coefficient vector ã_1,(k)⁺ is further normalized as indicated in line 7 of Table 1 to obtain a_1,(k+1)⁺.

TABLE 1
1: Input: a1,(0)+ = Initial random vector with unit norm
2: Input: μ > 0 = step size
3: Input: ε = Termination condition
4: k = 0
5: while ∥S11 · a1,(k)+∥22 > ε do
6: ã              1,(k)              + = a1,(k)+ − 2μS11H · S11 · a1,(k)+
7: a              1,(k+1)              + = ã1,(k)+/∥ã1,(k)+∥2
8: k = k + 1
9: end while

The first modal coefficient vector a_1,(0)⁺ at iteration k=0 may be the modal coefficient vector of any wavefront such that the modal coefficient vector has a unit norm. Alternatively, an additional wavefront may be selected and a first modal coefficient vector a_1,(0)⁺ be determined by normalizing the modal coefficient of the selected wavefront. The iteration repeats until the termination condition is reached, at which point the iterative refinement of the wavefront concludes. The final refined wavefront is taken as an approximation of the optimal wavefront a_opt. The termination condition occurs when the backscatter intensity 2 reaches or drops below a present level of ∈.

To physically implement the steepest descent algorithm, the gradient of the objective function may be rewritten as

S ₁₁ ^H ·S ₁₁ a _1,(k) ⁺ =F·S ₁₁ *·F·S ₁₁ ·a _1,(k) ⁺ =F·(S₁₁·(F·(S₁₁·a_1,(k)⁺)*))*  Equation 12

where F is defined as the matrix resulting from applying the operation flipud(*) to the identity matrix I (i.e., F=flipud(I)). The operation flipud(*) takes as its argument a vector or a matrix and flips the argument vector or matrix upside down such that the first row becomes the last row, the second row becomes the penultimate row, and so on. The above Equation 12 is physically realizable using the double phase conjugation method described below.

Double phase conjugation consists of (1) time reversing a wavefront, (2) transmitting the time reversed wavefront into the highly scattering random medium 100, and (3) time reversing the resulting backscatter wavefront. Where the starting wavefront used in the double phase conjugation method is the backscatter wavefront produced by transmitting the wavefront with modal coefficient vector a_1,(k)⁺ into the highly scattering random medium 100, double phase conjugation produces the objective function gradient as indicated in Equation 12. This process depends upon three matrix-vector operations: S₁₁·a₁⁺, F·(a₁⁻)*, and S₁₁^H·a₁⁻. Each of these matrix-vector operations may be performed mathematically or physically. In some embodiments, the example system 200 in FIG. 2 may be used to perform these operations.

FIG. 2 illustrates an example system for increasing transmission through the highly scattering random medium 100. The example system 200 may be used to iteratively refine a wavefront for transmission into the highly scattering random medium 100. In an embodiment of the system 200, a spatial light modulator (SLM) 202 is used to modulate an input wave 204 to produce a wavefront 206. Alternatively, a wavefront 206 may be produced directly or indirectly by any transmitter or modulator, using any technique known in the art. The wavefront 206 may then be focused on the highly scattering random medium 100 by passing through a lens 208. In some embodiments, the resulting backscatter wavefront 210 may pass back through the lens 208 and be redirected using a beam splitter 212. The backscatter wavefront 210 may then be time reversed as discussed below using a phase-conjugate mirror (PCM) 214. In some embodiments, the resulting time-reversed backscatter wavefront 216 may then be redirected using beam splitter 212, focused using lens 208, and transmitted into the highly scattering random medium 100 to produce another backscatter wavefront 210, which may pass back through the lens 208, be redirected using a beam splitter 212, and be time reversed using PCM 214 to produce another time-reversed backscatter wavefront 216. In some embodiments, the wavefront 216 may be redirected using a beam splitter 218 to reach a detector 220. The detector 220 may be a charge-coupled device or any other analog or digital transducer capable of detecting the wavefront 216. In some embodiments, the detector 220 may also receive a holographic reference wave 222. The holographic reference wave 222 may be used to determine the modal coefficient vector of the wavefront 216 using known digital holography techniques. Information corresponding to the wavefronts 216 and 222 detected at the detector 220 may be transmitted to a controller 224. The transmitted information may include the raw or processed measurements at the detector 220, or it may include one or more modal coefficient vectors or other identifying information associated with the wavefronts. In some embodiments, the controller 224 may be configured to reconstruct the wavefront 216 from the information received from the detector 220. In some embodiments, the controller 224 may perform computational matrix or vector operations to the modal coefficient vector of wavefront 216 to determine a new wavefront. The controller may then control the SLM 202 or other wave source to produce another wavefront 206 for transmission into the highly scattering random medium 100. Thus, the system 200 may be used to refine a wavefront by an iterative process to increase transmission through the highly scattering random medium 100. In some embodiments, the system 200 may also be used to focus a wavefront using an iterative process to increase field intensity at a point ρ.

Operation S₁₁a₁⁺ may be physically implemented by transmitting a wavefront 206 with modal coefficient vector a₁⁺ into the highly scattering random medium 100 and measuring (not shown) the resulting backscatter wavefront 210 with coefficient vector a₁⁻. The modal coefficient vector a₁⁻ may be determined from backscatter intensity measurements using known digital holography techniques implemented on the controller 224. In some embodiments, the modal coefficient vector a₁⁻ may not be determined; the wavefront 210 with modal coefficient vector a₁⁻ may instead be utilized without determining the vector a₁⁻, as below.

Operation F·(a₁⁻)* may be physically implemented by time reversing the wavefront 210 with coefficient vector a₁⁻. In some embodiments, the wavefront 210 may be time reversed by phase conjugation mirroring using PCM 214. Alternatively, the wavefront may be time reversed by determining the coefficient vector a₁⁻ and constructing a time reversed wavefront using the controller 224 by calculating the product of F and the complex conjugate of a₁⁻.

Operation S₁₁^H·a₁⁻ may be physically implemented by the double phase conjugation method by time reversing the wavefront 210 with modal coefficient a₁⁻, transmitting the time-reversed wavefront 216 into the highly scattering random medium 100, and time reversing the resulting backscatter wavefront 210. Where the highly scattering random medium 100 exhibits reciprocity, the reflection matrix S₁₁ must be such that the following relationship holds: S₁₁^H=F·S₁₁*·F. As a result, if transmitting a wave a into a highly scattering random medium 100 with reflection matrix S₁₁ produces a backscatter wave b, then transmitting a wave F·(a)* into a highly scattering random medium 100 with reflection matrix S₁₁^H produces a backscatter wave F·(b)*. From this relationship, the operation S₁₁^H·a₁⁻ may be expressed as F·S₁₁*, F·a₁⁻=F·(S₁₁·(F·(a₁⁻)*))*. The operation F·(a₁⁻)* corresponds to a time reversal of the wavefront 210 with modal coefficient vector a₁⁻, so the double phase conjugation method applied to wavefront 210 with modal coefficient vector a₁⁻ obtained by transmitting wavefront 206 into the highly scattering random medium 100 corresponds to the matrix-vector operation S₁₁^H·a₁⁻, each step of which may be implemented as described above.

Table 2 illustrates the steps of the kth iteration of the steepest descent method as physically implemented by applying the double phase conjugation method and the operations described above. In line 1, a first transmission wavefront 206 with modal coefficient vector a_1,(k)⁺ is transmitted into the highly scattering random medium 100, and the first backscatter wavefront 210 with modal coefficient vector a₁⁻ is obtained. Lines 2-4 correspond to the double phase conjugation method applied to the first backscatter wavefront 210 with modal coefficient vector a₁⁻. In line 2, the first backscatter wavefront 210 with modal coefficient vector a₁⁻ may be time reversed using PCM 214 to obtain the second transmission wavefront 216 with modal coefficient vector 4. In line 3, the second transmission wavefront 216 with modal coefficient vector a₁⁺ is transmitted into the highly scattering random medium 100, and a second backscatter wavefront 210 with modal coefficient vector a₁⁻ is obtained, which will not generally be the same as the first backscatter wavefront 210 with modal coefficient vector a₁⁻ obtained in line 1. In line 4, the second backscatter wavefront 210 with modal coefficient vector a₁⁻ may be time reversed using PCM 214 to obtain the updating wavefront 216 with modal coefficient vector 4, completing the double phase conjugation. In line 5, the wavefront a_1,(k)⁺ is refined in proportion to the step size μ and the modal coefficient vector a₁⁺ of updating wavefront 216 resulting from the double phase conjugation method to obtain a wave with modal coefficient vector ã₁⁺. The modal coefficient vector a₁⁺ may be determined using known digital holography techniques implemented on the detector 220, the holographic reference wave 222, and the controller 224. Finally, in line 6, the refined wavefront ã₁⁺ is normalized to obtain a_1,(k+1)⁺. The refinement and normalization in lines 5 and 6 may be performed using the controller 224. The iterative process in Table 2 terminates when the backscatter intensity of the first backscatter wavefront a₁⁻ falls below the present level of ∈, as indicated by ∥a₁⁻∥₂²<∈, which may also be determined using the controller 224.

TABLE 2
Vector Operation                 Physical Operation
1: α1− = S11· α1,(k)+            1                                                                    :                                                                                                                                                                                                              α                          _                                                                          1                          ,                                                      (                            k                            )                                                                          +                                                                                                        →                                              B                                                ackscatter                                                                                                                                α                        _                                            1                      -
2: α1+ = F · (α1−)*              2                                                                    :                                                                                                                                                                                                              α                          _                                                1                        -                                                                                                        →                      PCM                                                                                                          α                        _                                            1                      +
3: α1− = S12· α1+                3                                                                    :                                                                                                                                                                                                              α                          _                                                1                        +                                                                                                        →                      Backscatter                                                                                                          α                        _                                            1                      -
4: α1+ = F · (α1−)*              4                                                                    :                                                                                                                                                                                                              α                          _                                                1                        -                                                                                                        →                      PCM                                                                                                          α                        _                                            1                      +
5: {tilde over (α)}1+ = α1,(k)+ − 2μα1+ 5: {tilde over (α)}1+ = α1,(k)+ − 2μα1+
6: α1,(k+1)+ = α1+/∥α1+∥2        6                                                                    :                                                                                                                                                                                                                                          α                            _                                                    _                                                1                        +                                                                                                        →                      Normalization                                                                                                          α                        _                                                                    1                        ,                                                  (                                                      k                            +                            1                                                    )                                                                    +

The steepest descent method may also be implemented in embodiments comprising phase-only modulation of the wavefronts. As discussed above, phase-only modulation requires that all transmitted wavefronts have incident power ∥a₁⁺∥₂=1, or a₁⁺∈P₁^M. Table 3 illustrates the steepest descent method adjusted for phase-only modulation. In addition to the steps described above, the first backscatter wavefront 210 with modal coefficient vector a₁⁻ is used to determine a scalar a of the average magnitude of the modal coefficients of the first backscatter wavefront 210 with modal coefficient vector a₁⁻ as follows:

$\begin{array}{cc}\stackrel{\_}{a}=\frac{\sum _{n=-N}^N{\underset{\_}{a}}_{1,n}^{-}}{\sqrt{M}}& \mathrm{Equation}13\end{array}$

TABLE 3
Vector Operation                 Physical Operation
1: α1− = S11· α1,(k)+            1                                                                    :                                                                                                                                                                                                              α                          _                                                                          1                          ,                                                      (                            k                            )                                                                          +                                                                                                        →                      Backscatter                                                                                                          α                        _                                            1                      -
2                                                                    :                                                                                                                                                                                    α                        _                                                              =                                                                                            Σ                                                      n                            =                                                          -                              N                                                                                N                                                |                                                                              α                                                          1                              ,                              n                                                                                _                                                |                                                                    M 2                                                                    :                                                                                                                                                                                    α                        _                                                              =                                                                                            Σ                                                      n                            =                                                          -                              N                                                                                N                                                |                                                  α                                                      1                            ,                            n                                                    -                                                |                                                                    M
3: α1− ← p(/α1−)                 3: α1− ← p(/α1−)
4: α1+ = F · (α1−)*              4                                                                    :                                                                                                                                                                                                              α                          _                                                1                        -                                                                                                        →                      PCM                                                                                                          α                        _                                            1                      +
5: α1− = S11 · α1+               5                                                                    :                                                                                                                                                                                                              α                          _                                                1                        +                                                                                                        →                      Backsetter                                                                                                          α                        _                                            1                      -
6: α1+ = F · (α1−)*              6                                                                    :                                                                                                                                                                                                              α                          _                                                1                        -                                                                                                        →                      PCM                                                                                                          α                        _                                            1                      +
7: {tilde over (α)}1+ = α1,(k)+ − 2μ αα1+ 7: {tilde over (α)}1+ = α1,(k)+ − 2μ αα1+
8: α1,(k+1)+ = p(/{tilde over (α)}1+) 8: α1,(k+1)+ = p(/{tilde over (α)}1+)

This may be implemented by detecting the first backscatter wavefront 216 with detector 220, then computing the average magnitude scalar ā using the controller 224. Using this physically measurable scalar, the gradient of the objective function S₁₁^H·S₁₁·a_1,(k)⁺ is physically implementable in a method or system utilizing phase-only modulation as āS₁₁^H·p(∠(S₁₁·a_1,(k)⁺)). Once the average magnitude scalar ā has been determined in line 2, the phase-only modal coefficient vector of the first backscatter wavefront with modal coefficient vector p(∠a₁⁻) is determined in line 3. This may be implemented by the controller 224, which may control the SLM 202 to produce the phase-only first backscatter wavefront 206 with modal coefficient vector p(∠a₁⁻). The double phase conjugation method is then applied to the phase-only first backscatter wavefront 206 with modal coefficient vector p(∠a₁⁻) using beam splitters 226 and 228 to direct the wave 206 to the PCM 214. The updating wavefront 216 with modal coefficient vector a₁⁺ may be obtained by implementing the double phase conjugation method as above starting with the phase-only first backscatter wavefront 206. The coefficient vector a₁⁺ of updating wavefront 216 may be determined using the detector 220 and controller 224, as above. In line 7, the controller 224 refines the wavefront a_1,(k)⁺ in proportion to the step size μ, the average magnitude scalar ā, and the updating wavefront modal coefficient vector a₁⁺ to obtain a modal coefficient vector ã₁⁺. To satisfy the incident power constraint, the controller 224 uses only the phase of ã₁⁺ to produce the wavefront a_1,(k+1)⁺=p(∠ã₁⁺)∈P₁^M using the SLM 202 in line 8. The iterative process in Table 3 terminates when the backscatter intensity of the first backscatter wavefront a₁⁻ falls below the present level of ∈, as above.

The steepest descent method using phase-only modulation described above first determines an iteratively refined approximation ã₁⁺ of the optimal wavefront a_opt and then projects it onto the set of phase-only modulated wavefronts a_1,(k+1)⁺=p(∠ã₁⁺)∈P₁^M. Table 4 illustrates the phase-only gradient method, which updates only the phase of the incident wavefront θ_1,(k)⁺. to iteratively approximate the optimal phase vector θ_opt. In this method, the objective function of interest is the intensity of the backscatter of a phase-only modulated wavefront ∥S₁₁·p(θ)∥₂², the negative gradient of which is used to refine the phase of the incident wavefront θ_1,(k)⁺ in each iteration as follows:

$\begin{array}{cc}{\underset{\_}{\theta }}_{1,\left(k+1\right)}^{+}={\underset{\_}{\theta }}_{1,\left(k\right)}^{+}-\sqrt{M}\mu \frac{\partial {S_{11}·\underset{\_}{p}\left(\underset{\_}{\theta }\right)}_2^2}{\partial \underset{\_}{\theta }}{|}_{\underset{\_}{\theta }={\underset{\_}{\theta }}_{1,\left(k\right)}}.& \mathrm{Equation}14\end{array}$

This equation may be rewritten as follows:

θ _1,(k+1) ⁺=θ_1,(k)⁺−2√√{square root over (M)}μIm[diag{{circumflex over (p)}(−θ_1,(k)⁺)}·S₁₁^H·S₁₁·p(θ_1,(k)⁺)].  Equation 15

where Im(•) denotes the operator that returns the imaginary part of the argument and diag{p(−θ_1,(k)⁺)} denotes a diagonal matrix with entries p(θ_1,(k)⁺) along its diagonal. As above, the incident power constrain prevents the physical implementation of S₁₁^H·S₁₁·p(θ_1,(k)⁺), so āS₁₁^H·p(∠(S₁₁·p(θ_1,(k)⁺))) is substituted where ā is the average magnitude scalar defined above. Therefore, iterative refinement of the wavefront according to Equation 15 above can be rewritten in a physically implementable form as

θ _1,(k+1) ⁺=θ_1,(k)⁺−2√{square root over (M)}μ aIM[diag{_p(−θ_1,(k)⁺)}·S₁₁^H·p(∠S₁₁·p(θ_1,(k)⁺))],  Equation 16

TABLE 4
Vector Operation                 Physical Operation
1: α1− = S11· p(θ1,(k)+)         1                                                                    :                                                                                                                                                                                                              p                          _                                                                                                  (                                                                                    θ                              _                                                                                      1                              ,                                                              (                                k                                )                                                                                      +                                                    )                                                                                                                                →                      Backscatter                                                                                                          α                        _                                            1                      -
2                                                                    :                                                                                                                                                                                    α                        _                                                              =                                                                                            Σ                                                      n                            =                                                          -                              N                                                                                N                                                |                                                                              α                                                          1                              ,                              n                                                                                _                                                |                                                                    M 2                                                                    :                                                                                                                                                                                    α                        _                                                              =                                                                                            Σ                                                      n                            =                                                          -                              N                                                                                N                                                |                                                  α                                                      1                            ,                            n                                                    -                                                |                                                                    M
3: α1− ← p(/α1−)                 3: α1− ← p(/α1−)
4: α1+ = F · (α1−)*              4                                                                    :                                                                                                                                                                                    α                        1                        -                                                                                                        →                      PCM                                                                                                          α                        _                                            1                      +
5: α1− = S11 · α1+               5                                                                    :                                                                                                                                                                                                              α                          _                                                1                        +                                                                                                        →                      Backscatter                                                                                                          α                        _                                            1                      -
6: α1+ = F · (α1−)*              6                                                                    :                                                                                                                                                                                                              α                          _                                                1                        -                                                                                                        →                      PCM                                                                                                          α                        _                                            1                      +
7                                                                    :                                                                                                                                                                                                              θ                          _                                                                          1                          ,                                                      (                                                          k                              +                              1                                                        )                                                                          +                                                              =                                                                                            θ                          _                                                                          1                          ,                                                      (                            k                            )                                                                          +                                            -                                              2                                                                          M                                                                        μ                                                                          α                          _                                                                        1                                                                          m                                                                                [                                                          diag                                                                                                                              {                                                                                                            p                                      _                                                                                                                                                  (                                                                              -                                                                                                                              θ                                            _                                                                                                                                1                                            ,                                                                                          (                                              k                                              )                                                                                                                                +                                                                                                                    )                                                                                                        }                                                                ·                                                                                                      α                                    _                                                                    1                                  +                                                                                                                      ]

In line 1 of Table 4, a first transmission wavefront 206 with phase vector p(θ_1,(k)⁺) is transmitted into the highly scattering random medium 100, and the first backscatter wavefront 210 with modal coefficient vector a₁⁻ is obtained. Lines 2-6 determine the average magnitude scalar ā and implement the double phase conjugation method as above. In line 7, the controller 224 refines the phase of the incident wavefront θ_1,(k)⁺ according to Equation 16, using the average magnitude scalar ā and the results of the double phase conjugation process. The refined phase vector θ_{1, (k+1)}⁺ may be directly computed using the controller 224. The iterative process in Table 4 terminates when the backscatter intensity of the first backscatter wavefront a₁⁻ falls below the present level of ∈, as above.

The methods described above amplify the highly backscattering component of the wavefront 206 with modal coefficient vector a_1,(k)⁺ or p(θ_1,(k)⁺) in a manner analogous to time-reversal focusing, allowing the highly backscattering component to be iteratively removed from the wavefront. The refined wavefront will then produce less backscatter, which makes measurement of the backscatter modal coefficient vector increasingly difficult in subsequent iterations of the method. Furthermore, the step size μ must be set to a value that ensures convergence of the method. The conjugate gradient method described below maintains a high backscatter field intensity and does not require the selection of a step size.

Another method of determining the optimal wavefront a_opt involves a conjugate gradient algorithm as set forth in Table 5. The conjugate gradient method involves iteratively determining the backscatter minimizing wavefront

a _1,(k+1) ⁺ =a _1,(k) ⁺+μ_(k+1)d_(k),  Equation 17

where μ_(k+1) is a step size and d_(k) is a search direction. Specifically, the step size is determined as

μ_(k+1)=∥r_(k)∥₂²/∥S₁₁·d_(k)∥₂²  Equation 18

and the search direction is determined as

d _(k+1) =r _(k+1)+β_(k+1)d_(k)  Equation 19

where

β_(k+1)=∥r_(k+1)∥₂²/∥r_(k)∥₂²  Equation 20

and the residual vector r_(k+1) is determined as the physically realizable product

r _(k+1) =−S ₁₁ ^H ·S ₁₁ ·a _1,(k+1) ⁺.  Equation 21

Analogous to the steepest descent method above, the conjugate gradient method uses a termination condition such that iterative refinement terminates when ∥k_(k+1)∥₂<∈, which may be implemented at controller 224. The initial modal coefficient vector as a_a,(0)⁺ may be any unit vector, and d₍₀₎ and {circumflex over (r)}₍₀₎ may both be set to −S₁₁^H·S₁₁·a_1,(0)⁺ by application of the double phase conjugation method applied to the backscatter wavefront a₁⁻ obtained by transmitting −a_1,(0)⁺ into the highly scattering random medium 100.

Table 5 illustrates the vector and physical operations of each iteration of the conjugate gradient method. It should be noted that modal coefficient vector a_1,(k)⁺ is refined in each iteration at the controller 224, but the wavefront with modal coefficient vector a_1,(k)⁺ is not transmitted into the highly scattering random medium 100 until the termination condition has been met. This has the advantage of maintaining a higher backscatter field intensity throughout the iterative process, thereby reducing the role of noise and measurement error in the system. In line 1, a first transmission wavefront 206 with modal coefficient vector d_(k) is generated and transmitted into the highly scattering random medium 100, and the first backscatter wavefront 210 with modal coefficient vector d₁⁻ is obtained. Lines 2-4 correspond to the double phase conjugation method applied to the first backscatter wavefront d₁⁻. In line 2, the first backscatter wavefront 210 with modal coefficient vector d₁⁻ is time reversed using PCM 214 or another method to obtain the second transmission wavefront 216 with modal coefficient vector d₁⁺. In line 3, the second transmission wavefront 216 with modal coefficient vector d₁⁺ is transmitted into the highly scattering random medium 100, and a second backscatter wavefront 210 with modal coefficient vector di is obtained, which will not generally be the same as the first backscatter wavefront 210 with modal coefficient vector d₁⁻ obtained in line 1. In line 4, the second backscatter wavefront 210 with modal coefficient vector d₁⁻ is time reversed using PCM 214 or another method to obtain the updating wavefront 216 with modal coefficient vector d, completing the double phase conjugation. The modal coefficient vector {circumflex over (d)} of the updating wavefront 216 may be determined using the detector 220, holographic reference wave 222, and controller 224 using known digital holography techniques. In lines 5-8, the updated step size μ_(k+1), residual vector r_(k+1), factor β_(k+1), and search direction d_(k+1) are determined using the modal coefficient vector d. These updates may be directly or indirectly computed by standard matrix or vector operations using the controller 224. Finally, line 9 indicates the iterative refinement of the wavefront a_1,(k)⁺ using the updated step size μ_(k+1) and the search direction d_(k), and line 10 indicates the normalization of the refined wavefront a_1,(k+1)⁺. Both steps 9 and 10 may be implemented using controller 224. When the termination condition is achieved, iterative refinement terminates, and the refined wavefront 206 with modal coefficient vector a_1,(k+1)⁺ is transmitted into the highly scattering random medium 100.

TABLE 5
Vector Operation	Physical Operation
1:	d 1 − = S11 · d(k)	1:	${\underset{\_}{d}}_{\left(k\right)}\stackrel{\mathrm{Backscatter}}{\to }{\underset{\_}{d}}_1^{-}$
2:	d 1 + = F · (α1−)*	2:	${\underset{\_}{d}}_1^{-}\stackrel{\mathrm{PCM}}{\to }{\underset{\_}{d}}_1^{+}$
3:	d 1 − = S11 · d1+	3:	${\underset{\_}{d}}_1^{+}\stackrel{\mathrm{Backscatter}}{\to }{\underset{\_}{d}}_1^{-}$
4:	d = F · (d1−)*	4:	${\underset{\_}{d}}_1^{-}\stackrel{\mathrm{PCM}}{\to }\underset{\_}{d}$
5:	μ(k+1) = ∥r(k)∥22/(d(k)H · d)	5:	μ(k+1) = ∥r(k)∥22/(d(k)H · d)
6:	r (k+1) = r(k) − μ(k+1)d	6:	r (k+1) = r(k) − μ(k+1)d
7:	β(k+1) = ∥r(k+1)∥22/∥r(k)∥22	7:	β(k+1) = ∥r(k+1)∥22/∥r(k)∥22
8:	d (k+1) = r(k+1) + β(k+1)d(k)	8:	d (k+1) = r(k+1) + β(k+1)d(k)
9:	α 1,(k+1) + = α1,(k)+ + μ(k+1)d(k)	9:	α 1,(k+1) + = α1,(k)+ + μ(k+1)d(k)
10:	α 1,(k+1) + = α1,(k+1)+/∥α1,(k+1)+∥2	10:	α 1,(k+1) + = α1,(k+1)+/∥α1,(k+1)+∥2

Focusing Transmission Through Highly Scattering Random Media

The methods discussed above may also be applied to focus transmission through the highly scattering random medium 100 to a location ρ₀ in the halfspace where y>L. Since a₂⁻=0, the transmitted wave is a₂⁺=S₂₁·a₁⁺. Combining this with Equation 1, the field at location ρ₀ is

$\begin{array}{cc}{e}_2^{+}\left({\underset{\_}{\rho }}_{0}\right)=\underset{=:{\underset{\_}{f}\left({\underset{\_}{\rho }}_0\right)}^{″}}{\underset{}{\left[h_{-N}{}^{-j{\underset{\_}{k}}_N^{+}·{\rho }_{0}}\dots h_{-N}{}^{-j{\underset{\_}{k}}_N^{+}·{\underset{\_}{\rho }}_0}\right]}}·S_{21}·{\underset{\_}{a}}_1^{+}.& \mathrm{Equation}22\end{array}$

so the optimal focusing wavefront with modal coefficient vector a_foc may be determined as the power-constrained wavefront that maximizes the intensity of the field at location ρ₀ by transmission through the highly scattering random medium 100 with transmission matrix S₂₁ as follows:

$\begin{array}{cc}{\underset{\_}{a}}_{\mathrm{foc}}=\underset{{\underset{\_}{a}}_1^{+}}{\arg \max }\frac{{e_2^{+}\left({\underset{\_}{\rho }}_0\right)}_2^2}{{{\underset{\_}{a}}_1^{+}}_2^2}=\underset{{{\underset{\_}{a}}_1^{+}}_2=1}{\arg \max }{\underset{\underset{=:{\underset{\_}{c}\left({\underset{\_}{\rho }}_0\right)}^{″}}{}}{{\underset{\_}{f}}^H\left({\underset{\_}{\rho }}_0\right)·S_{21}}·a_1^{+}}_2^2,& \mathrm{Equation}23\end{array}$

the solution to which is

$\begin{array}{cc}{\underset{\_}{a}}_{\mathrm{foc}}=\frac{\underset{\_}{c}\left({\underset{\_}{\rho }}_0\right)}{{\underset{\_}{c}\left({\underset{\_}{\rho }}_0\right)}_2}=\frac{S_{21}^H·\underset{\_}{f}\left({\underset{\_}{\rho }}_0\right)}{{S_{21}^H·\underset{\_}{f}\left({\underset{\_}{\rho }}_0\right)}_2}.& \mathrm{Equation}24\end{array}$

The optimal focusing wavefront may be determined by time reversing the wavefront received by placing a wave source at location ρ₀. This has motivated work that determines the optimal focusing wavefront using O(M) measurements of the intensity at location ρ₀. The method described herein, however, uses the K eigen-wavefronts with transmission coefficients near unity (i.e., τ≈1) to determine the optimal focusing wavefront with modal coefficient vector a_foc using on the order of K measurements of the intensity at location ρ₀.

Expressing the optimal focusing wavefront in terms of the singular value decomposition of transmission matrix S₂₁ in Equation 24 above yields the expression

$\begin{array}{cc}{\underset{\_}{a}}_{\mathrm{foc}}\propto S_{21}^H·\underset{\_}{f}\left({\underset{\_}{\rho }}_{0}\right)=\sum _{i=1}^M{\sigma }_i\underset{\underset{=:w_i}{}}{\left({\underset{\_}{u}}_i^H·\underset{\_}{f}\left({\underset{\_}{\rho }}_0\right)\right)}{\underset{\_}{v}}_i=\sum _{i=1}^M{\sigma }_iw_i{\underset{\_}{v}}_i.& \mathrm{Equation}25\end{array}$

where σ_i is the singular value associated with the left and right singular vectors u_i and v_i, respectively. By convention, the singular values are arranged in descending orders so that σ₁≧σ₂≧ . . . ≧σ_M. The transmission coefficient for each right singular vector is τ(v_i)=σ_i², so only the K eigen-wavefronts with high transmission will significantly contribute to the field at location ρ₀. The remaining right singular vectors will have transmission coefficients close to zero. Since there are typically K>>M highly transmitting eigen-wavefronts, determining the optimal focusing wavefront with modal coefficient vector a_foc as a superposition of these K eigen-wavefronts of transmission matrix S₂₁ will produce an approximation of the optimal focusing wavefront using fewer measurements—O(K)>>O(M). As before, these K high-transmission eigen-wavefronts of transmission matrix S₂₁ correspond precisely to the K low-reflection eigen-wavefronts of reflection matrix S₁₁ in a lossless setting. Therefore, the optimal focusing wavefront with modal coefficient vector a_foc may be approximated using O(K) backscatter measurements to determine the backscattering-minimizing eigen-wavefronts of reflection matrix S₁₁ and O(K) intensity measurements at location ρ₀.

To determine the optimal focusing wavefront in an iterative manner, a generalized coordinate descent method for amplitude and phase optimization is used. This involves an M×N_B matrix B=[b₁ . . . b_NB] with orthonormal columns, such that B^H·B=I_NB. Thus, N_B denotes the number of orthonormal basis vectors. The modal coefficient vector a₁⁺ in Equation 24 can then be expanded in terms of the basis vectors given by the columns of B as

$\begin{array}{cc}{\underset{\_}{a}}_1^{+}=\sum _{i=1}^{N_B}p_i{}^{{\mathrm{j\phi }}_i}{\underset{\_}{b}}_i,& \mathrm{Equation}26\end{array}$

where p_l≧0 and φ_l∈[−π, π] are the unknown amplitudes and phases, respectively. The amplitudes p_l may be approximated by transmitting a wavefront 206 with modal coefficient vector a₁⁺=b_l into the highly scattering random medium 100 for each l=1, . . . , N_B, measuring the intensity of the field I_l at the location ρ₀ using a charge-coupled device or other detector (not shown), and setting p_l=√{square root over (I_l)} using the controller 224. The phases φ_l may similarly be approximated by setting each phase φ_l to any starting value, then sequentially finding the phase that optimizes the intensity of the field I_l at the location ρ₀ for each l=1, . . . , N_B. This may be accomplished using a simple line search to scan the intensity measured at location ρ₀ over a fixed set of discretized values of the phase φ_l, which may be implemented by the controller 224 controlling the SLM 202. Alternatively, other known or later developed algorithms to approximate each phase φ_l. Both of the processes of determining the amplitudes p_l and the phases φ_l may be accomplished using O(N_B) measurements. The number of measurements can be reduced to O(K) by using the K right singular eigenvectors v_i of reflection matrix S₁₁ that minimize backscattering as the N_B columns of matrix B.

The K right singular eigenvectors v_i of reflection matrix S₁₁ that minimize backscattering correspond to the right singular eigenvectors v_i associated with the K smallest singular values σ_i of reflection matrix S₁₁. The K smallest singular values σ_i and right singular vectors v_i may be determined using the Lanczos algorithm. The Lanczos algorithm iteratively produces a tridiagonal matrix H, the eigenvalues and eigenvectors of which approximate the eigenvalues and eigenvectors of S₁₁·S₁₁. Table 6 illustrates an iteration of a physically implementable method corresponding to the Lanczos Algorithm that utilizes double phase conjugation, as described above. In the first iteration where k=1, the modal coefficient vector of the iterative wavefront q_(k) may be set to any unit norm vector and the coefficient s₍₀₎=0. In line 1, the iterate wavefront 206 with modal coefficient vector q_(k) is transmitted into the highly scattering random medium 100, and the backscatter wavefront 210 with modal coefficient vector q₁⁻ is obtained. Lines 2-4 correspond to the double phase conjugation method applied to the first backscatter wavefront 210 with modal coefficient vector q₁⁻ to obtain the wavefront 216 with modal coefficient vector v, which may be determined using controller 224 from measurements at detector 220. In line 5, the kth diagonal entry H_k,k in tridiagonal matrix H is determined as the product of the complex conjugate transpose of the modal coefficient vector q_(k) of the iterate wavefront 206 and modal coefficient vector v. Line 6 updates the modal coefficient vector v using the product of the diagonal entry H_k,k and the modal coefficient vector {circumflex over (q)}_(k) of the iterate wavefront 206 and the product of the previous iteration coefficient s_(k-1) and the modal coefficient vector q_(k−1) of the previous iterative wavefront 206. In line 7, the entries in the diagonals of matrix H above and below the main diagonal are determined as the norm of modal coefficient vector v, as is coefficient s_(k). Finally, in line 8, the modal coefficient vector q_(k+1) of the next iterative wavefront 206 is determined as the quotient of updated modal coefficient vector v and coefficient s_(k). The matrix and vector operations in lines 5-8 may be directly or indirectly computed using the standard matrix or vector operations implemented using the controller 224.

TABLE 6
Vector Operation	Physical Operation
1:	q 1 − = S11 · q(k)	1:	${\underset{\_}{q}}_{\left(k\right)}\stackrel{\mathrm{Backscatter}}{\to }{\underset{\_}{q}}_1^{-}$
2:	q 1 + = F · (q1−)*	2:	${\underset{\_}{q}}_1^{-}\stackrel{\mathrm{PCM}}{\to }{\underset{\_}{q}}_1^{+}$
3:	q 1 − = S11 · q1+	3:	${\underset{\_}{q}}_1^{+}\stackrel{\mathrm{Backscatter}}{\to }{\underset{\_}{q}}_1^{-}$
4:	v = F · (q1−)*	4:	${\underset{\_}{q}}_1^{-}\stackrel{\mathrm{PCM}}{\to }\underset{\_}{v}$
5:	Hk,k = q(k)H · v	5:	Hk,k = q(k)H · v
6:	v = v − Hk,kq(k) − s(k−1)q(k−1)	6:	v = v − Hk,kq(k) − s(k−1)q(k−1)
7:	Hk+1,k = Hk,k+1 = s(k) = ∥v∥2	7:	Hk+1,k = Hk,k+1 = s(k) = ∥v∥2
8:	q (k+1) = v/s(k)	8:	q (k+1) = v/s(k)

From the Lanczos algorithm as implemented above, the matrix B may be determined using controller 224 as

B=Q·U  Equation 27

where Q=[q₍₁₎ . . . {circumflex over (q)}_(NB))] is the matrix with the modal coefficient vectors of the iterative wavefronts q_(k) as its columns and U=[u₍₁₎ . . . u_(NB)] is the matrix with the N_B eigenvectors of matrix H associated with the N_B smallest eigenvalues of matrix H as its columns. The convergence theory of the Lanczos algorithm predicts that the eigenvector estimates will rapidly converge to the K eigenvectors of S₁₁^H·S₁₁ associated with the eigen-wavefronts of S₁₁ with the smallest reflection coefficients; hence, setting N_B to be O(K) will suffice. The number of highly transmitting right eigenvectors K may be estimated by counting the number of eigenvalues of matrix H that are below a present level of ∈. Therefore, the matrix B will be of dimensions M×K, and the process described above for determining the amplitudes p_l and the phases φ_l of the optimal focusing wavefront with modal coefficient vector a_foc may be accomplished using O(K) measurements.

Some or all calculations performed by the controller 224 in the embodiments described above (e.g., calculations for determining the modal coefficient vector associated with a wavefront or matrix-vector operations) may be performed by a computer such as a personal computer, laptop computer, server, mainframe, or specialized computing device or system. The calculations performed by the controller 224 may also be performed by a program or programs implemented on one or more computers. In some embodiments, some or all of the calculations may be performed by more than one computer communicatively connected through a network.

FIG. 3 illustrates an exemplary block diagram of a computer system 300 on which an exemplary method for increasing transmission through a highly scattering random medium may operate in accordance with the described embodiments. The computer system 300 of FIG. 3 includes a computing device in the form of a computer 310. Components of the computer 310 may include, but are not limited to, a processing unit 320, a system memory 330, and a system bus 321 that couples various system components including the system memory to the processing unit 320. The system bus 321 may be any of several types of bus structures including a memory bus or memory controller, a peripheral bus, and a local bus using any of a variety of bus architectures. By way of example, and not limitation, such architectures include the Industry Standard Architecture (ISA) bus, Micro Channel Architecture (MCA) bus, Enhanced ISA (EISA) bus, Video Electronics Standards Association (VESA) local bus, and Peripheral Component Interconnect (PCI) bus (also known as Mezzanine bus).

Computer 310 typically includes a variety of computer readable media. Computer readable media can be any available media that can be accessed by computer 310 and includes both volatile and nonvolatile media, and both removable and non-removable media. By way of example, and not limitation, computer readable media may comprise computer storage media and communication media. Computer storage media includes volatile and nonvolatile, removable and non-removable media implemented in any method or technology for storage of information such as computer readable instructions, data structures, program modules or other data. Computer storage media includes, but is not limited to, RAM, ROM, EEPROM, FLASH memory or other memory technology, CD-ROM, digital versatile disks (DVD) or other optical disk storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other medium which can be used to store the desired information and which can accessed by computer 310. Communication media typically embodies computer readable instructions, data structures, program modules or other data in a modulated data signal such as a carrier wave or other transport mechanism and includes any information delivery media. The term “modulated data signal” means a signal that has one or more of its characteristics set or changed in such a manner as to encode information in the signal. By way of example, and not limitation, communication media includes wired media such as a wired network or direct-wired connection, and wireless media such as acoustic, radio frequency (RF), infrared and other wireless media. Combinations of any of the above are also included within the scope of computer readable media.

The system memory 330 includes computer storage media in the form of volatile and/or nonvolatile memory such as read only memory (ROM) 331 and random access memory (RAM) 332. A basic input/output system 333 (BIOS), containing the basic routines that help to transfer information between elements within computer 310, such as during start-up, is typically stored in ROM 331. RAM 332 typically contains data and/or program modules that are immediately accessible to and/or presently being operated on by processing unit 320. By way of example, and not limitation, FIG. 3 illustrates operating system 334, application programs 335, other program modules 336, and program data 337.

The computer 310 may also include other removable or non-removable, volatile or nonvolatile computer storage media. By way of example only, FIG. 3 illustrates a hard disk drive 341 that reads from or writes to non-removable, nonvolatile magnetic media, a magnetic disk drive 351 that reads from or writes to a removable, nonvolatile magnetic disk 352, and an optical disk drive 355 that reads from or writes to a removable, nonvolatile optical disk 356 such as a CD ROM or other optical media. Other removable or non-removable, volatile or nonvolatile computer storage media that can be used in the exemplary operating environment include, but are not limited to, magnetic tape cassettes, flash memory cards, digital versatile disks, digital video tape, solid state RAM, solid state ROM, and the like. The hard disk drive 341 is typically connected to the system bus 321 through a non-removable memory interface such as interface 340, and magnetic disk drive 351 and optical disk drive 355 are typically connected to the system bus 321 by a removable memory interface, such as interface 350.

The drives and their associated computer storage media discussed above and illustrated in FIG. 3 provide storage of computer readable instructions, data structures, program modules and other data for the computer 310. In FIG. 3, for example, hard disk drive 341 is illustrated as storing operating system 344, application programs 345, other program modules 346, and program data 347. Note that these components can either be the same as or different from operating system 334, application programs 335, other program modules 336, and program data 337. Operating system 344, application programs 345, other program modules 346, and program data 347 are given different numbers here to illustrate that, at a minimum, they are different copies. A user may enter commands and information into the computer 310 through input devices such as a keyboard 362 and cursor control device 361, commonly referred to as a mouse, trackball or touch pad. A monitor 391 or other type of display device is also connected to the system bus 321 via an interface, such as a graphics controller 390. In addition to the monitor, computers may also include other peripheral output devices such as printer 396, which may be connected through an output peripheral interface 395.

The computer 310 may operate in a networked environment using logical connections to one or more remote computers, such as a remote computer 380. The remote computer 380 may be a personal computer, a server, a router, a network PC, a peer device or other common network node, and typically includes many or all of the elements described above relative to the computer 310, although only a memory storage device 381 has been illustrated in FIG. 3. The logical connections depicted in FIG. 3 include a local area network (LAN) 371 and a wide area network (WAN) 373, but may also include other networks.

When used in a LAN networking environment, the computer 310 is connected to the LAN 371 through a network interface or adapter 370. When used in a WAN networking environment, the computer 310 typically includes a modem 372 or other means for establishing communications over the WAN 373, such as the Internet. The modem 372, which may be internal or external, may be connected to the system bus 321 via the input interface 360, or other appropriate mechanism. In a networked environment, program modules depicted relative to the computer 310, or portions thereof, may be stored in the remote memory storage device 381. By way of example, and not limitation, FIG. 3 illustrates remote application programs 385 as residing on memory device 381.

The communications connections 370, 372 allow the device to communicate with other devices. The communications connections 370, 372 are an example of communication media. The communication media typically embodies computer readable instructions, data structures, program modules or other data in a modulated data signal such as a carrier wave or other transport mechanism and includes any information delivery media. A “modulated data signal” may be a signal that has one or more of its characteristics set or changed in such a manner as to encode information in the signal. By way of example, and not limitation, communication media includes wired media such as a wired network or direct-wired connection, and wireless media such as acoustic, RF, infrared and other wireless media. Computer readable media may include both storage media and communication media.

Throughout this specification, plural instances may implement components, operations, or structures described as a single instance. Although individual operations of one or more methods are illustrated and described as separate operations, one or more of the individual operations may be performed concurrently, and may, where appropriate, be performed in an order other than the order illustrated. Structures and functionality presented as separate components in example configurations may be implemented as a combined structure or component. Similarly, structures and functionality presented as a single component may be implemented as separate components. These and other variations, modifications, additions, and improvements fall within the scope of the subject matter herein.

The various operations of example methods described herein may be performed, at least partially, by one or more processors that are temporarily configured (e.g., by software) or permanently configured to perform the relevant operations. Whether temporarily or permanently configured, such processors may constitute processor-implemented modules that operate to perform one or more operations or functions. The modules referred to herein may, in some example embodiments, comprise processor-implemented modules.

Similarly, the methods or routines described herein may be at least partially processor-implemented. For example, at least some of the operations of a method may be performed by one or more processors or processor-implemented hardware modules. The performance of certain of the operations may be distributed among the one or more processors, not only residing within a single machine, but deployed across a number of machines. In some example embodiments, the processor or processors may be located in a single location, while in other embodiments the processors may be distributed across a number of locations.

Unless specifically stated otherwise, discussions herein using words such as “processing,” “computing,” “calculating,” “determining,” “presenting,” “displaying,” or the like may refer to actions or processes of a machine (e.g., a computer) that manipulates or transforms data represented as physical (e.g., electronic, magnetic, or optical) quantities within one or more memories (e.g., volatile memory, non-volatile memory, or a combination thereof), registers, or other machine components that receive, store, transmit, or display information.

As used herein any reference to “one embodiment” or “an embodiment” means that a particular element, feature, structure, or characteristic described in connection with the embodiment is included in at least one embodiment. The appearances of the phrase “in one embodiment” in various places in the specification are not necessarily all referring to the same embodiment.

Some embodiments may be described using the expression “coupled” and “connected” along with their derivatives. For example, some embodiments may be described using the term “coupled” to indicate that two or more elements are in direct physical or electrical contact. The term “coupled,” however, may also mean that two or more elements are not in direct contact with each other, but yet still co-operate or interact with each other. The embodiments are not limited in this context.

As used herein, the terms “comprises,” “comprising,” “includes,” “including,” “has,” “having” or any other variation thereof, are intended to cover a non-exclusive inclusion. For example, a process, method, article, or apparatus that comprises a list of elements is not necessarily limited to only those elements but may include other elements not expressly listed or inherent to such process, method, article, or apparatus. Further, unless expressly stated to the contrary, “or” refers to an inclusive or and not to an exclusive or. For example, a condition A or B is satisfied by any one of the following: A is true (or present) and B is false (or not present), A is false (or not present) and B is true (or present), and both A and B are true (or present).

In addition, use of the “a” or “an” are employed to describe elements and components of the embodiments herein. This is done merely for convenience and to give a general sense of the description. This description, and the claims that follow, should be read to include one or at least one and the singular also includes the plural unless it is obvious that it is meant otherwise.

This detailed description is to be construed as exemplary only and does not describe every possible embodiment, as describing every possible embodiment would be impractical, if not impossible. One could implement numerous alternate embodiments, using either current technology or technology developed after the filing date of this application.

Claims

1. A computer-implemented method for transmission of a wavefront through a medium, the computer-implemented method comprising: refining, with one or more processors, the wavefront in each of one or more iterations, including: transmitting, from a transmitter, a test wavefront into the medium;receiving, at a detector, a measurement of the back-scatter from the test wavefront entering the medium;determining, with one or more processors, whether the test wavefront achieves a numerically prespecified level of transmission through the medium based upon the received back-scatter;refining, with one or more processors, the test wavefront for transmitting into the medium in the next iteration based upon the received back-scatter;terminating the iterations when the numerically prespecified level of transmission through the medium is achieved;modifying, with one or more processors, the wavefront based upon the final test wavefront; andtransmitting, from the transmitter, the modified wavefront through the medium.

2. The computer-implemented method of claim 1, wherein refining the test wavefront comprises one of the following: modifying the amplitude and phase of the test wavefront or modifying only the phase of the test wavefront.

3. The computer-implemented method of claim 1, wherein transmitting the test wavefront into the medium and receiving the measurement of the back-scatter from the test wavefront entering the medium in each of the one or more iterations comprises: transmitting a first test wavefront into the medium;receiving a measurement of a first backscatter wavefront from the first test wavefront entering the medium;time-reversing the first backscatter wavefront to obtain a second test wavefront;transmitting the second test wavefront into the medium;receiving a measurement of a second backscatter wavefront from the second test wavefront entering the medium; andtime-reversing the second backscatter wavefront to obtain a third test wavefront.

4. The computer-implemented method of claim 3, wherein refining the test wavefront for transmitting into the medium in the next iteration based upon the received back-scatter comprises adjusting the test wavefront in proportion to the third test wavefront and normalizing the adjusted test wavefront to have the same power as the wavefront.

5. The computer-implemented method of claim 3, wherein: time-reversing the first backscatter wavefront to obtain a second test wavefront further comprises determining an average magnitude of the modal coefficients of the first backscatter wavefront, determining a phase-modulated first backscatter wavefront with the same modal coefficients as the first backscatter wavefront, and time-reversing the phase-modulated first backscatter wavefront; andrefining the test wavefront for transmitting into the medium in the next iteration based upon the received back-scatter comprises adjusting the test wavefront in proportion to the product of the third test wavefront and the average magnitude of the modal coefficients of the first backscatter wavefront and determining a phase-modulated test wavefront with the same modal coefficients as the adjusted test wavefront.

6. The computer-implemented method of claim 3, wherein: the first backscatter wavefront in each of the one or more iterations is distinct from the wavefront; andrefining the test wavefront for transmitting into the medium in the next iteration based upon the received back-scatter comprises determining a residual vector from the third test wavefront.

7. The computer-implemented method of claim 3, wherein: time-reversing the first backscatter wavefront to obtain a second test wavefront further comprises determining an average magnitude of the modal coefficients of the first backscatter wavefront, determining a phase-modulated first backscatter wavefront with the same modal coefficients as the first backscatter wavefront, and time-reversing the phase-modulated first backscatter wavefront; andrefining the test wavefront for transmitting into the medium in the next iteration based upon the received back-scatter comprises adjusting the phase of the test wavefront in proportion to the average magnitude of the modal coefficients of the first backscatter wavefront.

8. The computer-implemented method of claim 3, wherein the time-reversing the first backscatter wavefront and the time-reversing the second backscatter wavefront are performed using a phase conjugate mirror.

9. The computer-implemented method of claim 3, wherein: each of the one or more iterations further comprises determining one or more entries in a matrix H based upon the first test wavefront and the third test wavefront;the computer-implemented method further comprises determining a matrix Q the columns of which are the modal coefficient vectors of the test wavefronts of each of the one or more iterations; andmodifying the wavefront based upon the final test wavefront further comprises determining an optimal focusing wavefront based upon the matrix Q and the eigenvectors of the matrix H.

10. A system for transmission of a wavefront through a medium, the system comprising: one or more processors;one or more spatial light modulators configured to modulate the wavefront and transmit the modulated wavefront into the medium;one or more detectors configured to measure the backscatter field resulting from the transmission of the wavefront into the medium; anda program memory storing non-transitory executable instructions that when executed by the one or more processors cause the system to: refine the wavefront in each of one or more iterations, including: causing the one or more spatial light modulators to transmit a test wavefront into the medium;measuring, with the one or more detectors, the back-scatter from the test wavefront entering the medium;determining, with the one or more processors, whether the test wavefront achieves a numerically prespecified level of transmission through the medium based upon the received back-scatter;refining, with the one or more processors, the test wavefront for transmitting into the medium in the next iteration based upon the received back-scatter;terminate the iterations when the numerically prespecified level of transmission through the medium is achieved;modify, with the spatial light modulator, the wavefront based upon the final test wavefront; andtransmit, with the spatial light modulator, the modified wavefront through the medium.

11. The system of claim 10, further comprising executable instructions that when executed by the processor cause the system to refine the wavefront in each of the one or more iterations by causing the one or more spatial light modulators to transmit a test wavefront into the medium and measuring, with the one or more detectors, the back-scatter from the test wavefront entering the medium further comprises: transmitting, with the spatial light modulator, a first test wavefront into the medium;receiving a measurement of a first backscatter wavefront from the first test wavefront entering the medium;time-reversing the first backscatter wavefront to obtain a second test wavefront;transmitting the second test wavefront into the medium;receiving a measurement of a second backscatter wavefront from the second test wavefront entering the medium; andtime-reversing the second backscatter wavefront to obtain a third test wavefront.

12. The system of claim 10, wherein refining, with the one or more processors, the test wavefront for transmitting into the medium in the next iteration based upon the received back-scatter further comprises adjusting, with the one or more processors, the test wavefront in proportion to the third test wavefront and normalizing, with the one or more processors, the adjusted test wavefront to have the same power as the wavefront.

13. The system of claim 11, wherein: the first backscatter wavefront in each of the one or more iterations is distinct from the wavefront; andrefining, with the one or more processors, the test wavefront for transmitting into the medium in the next iteration based upon the received back-scatter comprises determining, with the one or more processors, a residual vector from the third test wavefront.

14. The system of claim 11, further comprising one or more phase conjugate mirrors and wherein: time-reversing the first backscatter wavefront to obtain a second test wavefront further comprises time-reversing, with the one or more phase conjugate mirrors, the first backscatter wavefront to obtain a second test wavefront; andtime-reversing the second backscatter wavefront to obtain a third test wavefront further comprises time-reversing, with the one or more phase conjugate mirrors, the second backscatter wavefront to obtain a third test wavefront.

15. A tangible, non-transitory computer-readable medium storing instructions for optimizing transmission of a wavefront through a medium that, when executed by one or more processors of a computer system, cause the computer system to: refine the wavefront in each of one or more iterations, comprising: causing a test wavefront to be transmitted into the medium;receiving a measurement of the back-scatter from the test wavefront entering the medium;determining whether the test wavefront achieves a numerically prespecified level of transmission through the medium based upon the received back-scatter;refining the test wavefront for transmitting into the medium in the next iteration based upon the received back-scatter;terminate the iterations when the numerically prespecified level of transmission through the medium is achieved;cause the wavefront to be modified based upon the final test wavefront; andcause the modified wavefront to be transmitted through the medium.

16. The tangible, non-transitory computer-readable medium of claim 15, wherein the executable instructions that when executed by the one or more processors cause the computer system to refine the wavefront in each of one or more iterations by causing a test wavefront to be transmitted into the medium and receiving a measurement of the back-scatter from the test wavefront entering the medium further comprise executable instructions that when executed by the one or more processors cause the computer system to refine the wavefront in each of one or more iterations by: causing a first test wavefront to be transmitted into the medium;receiving a measurement of a first backscatter wavefront from the first test wavefront entering the medium;time-reversing the first backscatter wavefront to obtain a second test wavefront;causing the second test wavefront to be transmitted into the medium;receiving a measurement of a second backscatter wavefront from the second test wavefront entering the medium; andtime-reversing the second backscatter wavefront to obtain a third test wavefront.

17. The tangible, non-transitory computer-readable medium of claim 16, wherein: refining the test wavefront for transmitting into the medium in the next iteration based upon the received back-scatter comprises adjusting the test wavefront in proportion to the third test wavefront and normalizing the adjusted test wavefront to have the same power as the wavefront; andthe time-reversing the first backscatter wavefront and the time-reversing the second backscatter wavefront are performed using a phase conjugate mirror.

18. The tangible, non-transitory computer-readable medium of claim 16, wherein: time-reversing the first backscatter wavefront to obtain a second test wavefront further comprises determining an average magnitude of the modal coefficients of the first backscatter wavefront, determining a phase-modulated first backscatter wavefront with the same modal coefficients as the first backscatter wavefront, and time-reversing the phase-modulated first backscatter wavefront; andrefining the test wavefront for transmitting into the medium in the next iteration based upon the received back-scatter comprises adjusting the test wavefront in proportion to the product of the third test wavefront and the average magnitude of the modal coefficients of the first backscatter wavefront and determining a phase-modulated test wavefront with the same modal coefficients as the adjusted test wavefront.

19. The tangible, non-transitory computer-readable medium of claim 16, wherein: time-reversing the first backscatter wavefront to obtain a second test wavefront further comprises determining an average magnitude of the modal coefficients of the first backscatter wavefront, determining a phase-modulated first backscatter wavefront with the same modal coefficients as the first backscatter wavefront, and time-reversing the phase-modulated first backscatter wavefront; andrefining the test wavefront for transmitting into the medium in the next iteration based upon the received back-scatter comprises adjusting the phase of the test wavefront in proportion to the average magnitude of the modal coefficients of the first backscatter wavefront.

20. The tangible, non-transitory computer-readable medium of 16, wherein the executable instructions that when executed by the one or more processors cause the computer system to refine the wavefront in each of one or more iterations further comprise executable instructions that when executed by the one or more processors cause the computer system to determine one or more entries in a matrix H based upon the first test wavefront and the third test wavefront, wherein the executable instructions that when executed by the one or more processors cause the computer system to determine a matrix Q the columns of which are the modal coefficient vectors of the test wavefronts of each of the one or more iterations, and wherein the executable instructions that when executed by the one or more processors cause the computer system to cause the wavefront to be modified based upon the final test wavefront further comprise executable instructions that when executed by the one or more processors cause the computer system to determine an optimal focusing wavefront based upon the matrix Q and the eigenvectors of the matrix H and modify the wavefront based upon the optimal focusing wavefront.