Receiver estimation engine for a chaotic system

[0002] The present invention relates generally to chaotic communications systems and more particularly to a chaotic communication system utilizing a chaotic receiver estimation engine. This estimation engine both synchronizes and recovers data by mapping probability calculation results onto the chaotic dynamics via a strange attractor geometrical approximation. For synchronization, the system does not require either a stable/unstable subsystem separation or a chaotic system inversion. The techniques employed can be implemented with any chaotic system for which a suitable geometric model of the attractor can be found.

[0003] Chaos is an area of science and mathematics that is characterized by processes that are nonlinear, such as equations that have a squared or higher order term. Chaotic processes are iterative, in that they perform the same operations over and over. By taking the result of the process (equations) and performing the process again on the result, a new result is generated. Chaos is often called deterministically random motion. Chaotic processes are deterministic because they can be described by equations and because knowledge of the equations and a set of initial values allows all future values to be determined. Chaotic processes also appear random since a sequence of numbers generated by chaotic equations has the appearance of randomness. One unique aspect of chaos compared to an arbitrarily chosen aperiodic nonlinear process is that the chaotic process can be iterated an infinite number of times, with a result that continues to exist in the same range of values or the same region of space.

[0004] A chaotic system exists in a region of space called phase space. Points within the phase space fly away from each other with iterations of the chaotic process. Their trajetories are stretched apart but their trajectories are then folded back onto themselves into other local parts of the phase space but still occupy a confined region of phase space A geometrical shape, called the strange attractor, results from this stretching and folding process. One type of strange attractor for a chaotic process called the Rossler system is depicted in FIG. 1, and illustrates the stretching and folding process. The chaotic attracter exists in perpetuity in the region of phase space for which the chaotic system is stable.

[0005] In an unstable system, two points in phase space that are initially close together become separated by the stretching and the folding, which causes the points to be placed correspondingly farther apart in the phase space than they originally were. Repeating the process (iterating) accentuates the situation. The points not only diverge from each other, but both trajectories move toward infinity and, therefore, away from the stable region of perpetual chaotic existence.

[0006] The sequences of points that result from two closely spaced initial conditions become very different very quickly in chaotic processes. The Henon chaotic system, for example, has been shown to start with two initial conditions differring by one (1) digit in the 15^thdecimal place. The result was that within 70 iterates the trajectories were so different that subtracting them resulted in a signal that was as large as the original trajectories' signals themselves. Therefore, the stretching and folding causes the chaotic process to exhibit sensitive dependence on initial conditions. A receiver without perfect information about a transmitter (a “key”) will be unable to lock to the transmitter and recover the message. Even if a very close lock were achieved at one point in time, it is lost extremely quickly because the time sequences in the transmitter and receiver diverge from each other within a few iterations of the chaotic process.

[0007] Chaotic nonlinear dynamics may be utilized in telecommunications systems. There is interest in utilizing and exploiting the nonlinearities to realize secure communications, while achieving reductions in complexity, size, cost, and power requirements over the current communications techniques. Chaotic processes are inherently spread in frequency, secure in that they possess a low probability of detection and a low probability of intercept, and are immune to most of the conventional detection, intercept, and disruption methods used against current ‘secure’ communications systems based on linear pseudorandom noise sequence generators. Chaotic time sequences theoretically never repeat, making them important for such applications as cryptographic methods and direct sequence spread spectrum spreading codes. In addition, chaotic behavior has been found to occur naturally in semiconductors, feedback circuits, lasers, and devices operating in compression outside their linear region.

[0008] Unfortunately, many of these characteristics also complicate the task of recovering the message in the chaotic transmission. A fundamental problem permeating chaotic communications research is the need for synchronization of chaotic systems and/or probabilistic estimation of chaotic state, without which there can be no transfer of information. Without synchronization, meaning alignment of a local receiver chaotic signal or sequence with that of the transmitter, the characteristic of sensitive dependence on initial conditions causes free-running chaotic oscillators or maps to quickly diverge from each other, preventing the transfer of information. Probability-based estimates of chaotic state are typically made via correlation or autocorrelation calculations, which are time consuming. The number of iterates per bit can be reasonable in conjunction with synchronization, but values of several thousand iterates per bit are typically seen without synchronization.

[0009] A general communications system employing chaos usually consists of a message m(t) injected into a chaotic transmitter, resulting in a chaotic, encoded transmit signal y(t). This signal is altered as it passes through the channel, becoming the received signal r(t). The receiver implements some mix of chaotic and/or statistical methods to generate an estimate of the message m_e(t).

[0010] The recovery of a message from transmitted data depends on the ability of the receiver to perform either asynchronous or synchronous detection. There are three fundamental approaches to chaotic synchronization delineated in the literature of the prior art: (1) decomposition into subsystems, (2) inverse system approach, and (3) linear feedback.

[0011] If the approach is a decomposition into subsystems, the chaotic system is divided into two or more parts: an unstable subsystem usually containing a nonlinearity and one or more stable subsystems. The stable subsystem may contain a nonlinearity if its Lyapunov exponents remain negative. The Lyapunov exponents of all subsystems must be computed to evaluate the stability or instability of each one, where negative Lyapunov exponents indicate stability and positive Lyapunov exponents indicate instability. The arbitrary division into stable and unstable systems is accomplished by trial and error until successful. The driver or master system is the complete chaotic system, and the driven or slave system(s) consists of the stable subsystem(s). In general, synchronization will depend on circuit component choices and initial conditions.

[0012] The problem with using the decomposition approach is threefold. First, the decomposition into subsystems is arbitrary, with success being defined via the results of Lyapunov exponent calculations. There is presently no known method by which to synthesize successful subsystem decomposition to satisfy a set of requirements. Second, the evaluation of Lyapunov exponents can become extremely involved, complicated, and messy. The time and effort spent evaluating a trial, arbitrary decomposition can be extensive, with no guarantee of success. Third, there can be several possible situations for the subsystem decompositions for a given chaotic process. A chaotic process could have single acceptable decomposition, multiple acceptable decompositions, or no acceptable decompositions. Since the method is by trial-and-error and without guarantee of success, a general approach that will work for an arbitrary chaotic process is highly desirable.

[0013] In the inverse system approach, the receiver is designed to invert the chaotic function of the transmitter. For appropriate choices of transmitter systems and initial conditions, inverse systems can be designed into a receiver that results in perfect recovery of a noiseless transmitted message. However, one major drawback of using the inverse system approach is that signal noise results in erroneous message recovery. In addition, the invertible system must be based on a specially designed chaotic process since generic chaotic signals are not invertible. The collection of specially designed invertible chaotic processes is very small. This seriously limits the repertoire of receiver designs as well as the library of applicable chaotic process choices. Hostile countermeasures designers would have a relatively easy time analyzing and compromising the small set of systems that could be developed using these invertible chaotic functions.

[0014] In the linear feedback system approach, a master-slave setup is used in which the output of the slave system is fed back and subtracted from the received chaotic time sequence to generate an error function. The error function drives the slave system in a manner reminiscent of a phase-lock loop. Synchronization depends on the initial conditions of the master and slave systems. One drawback of linear feedback systems is that without a set of constraints on the initial conditions, synchronization of the master and slave systems cannot be guaranteed. In addition, the slave is a stable subsystem of the chaotic process, and is subject to the same drawbacks and limitations enumerated for the subsystem decomposition method above.

[0015] In addition to the three fundamental approaches to chaotic synchronization, there are four basic methods of encoding digital data into chaotic transmitted sequences: masking, orbit control, chaos shift key, and orthogonal autocorrelation sequence. They employ several different chaotic systems, with the most commonly mentioned including the Henon, Rossler, Lorenz, and Double Scroll (Chua's circuit). Most of these methods have variations in their implementations. A common problem among them is the issue of transmitter-receiver synchronization.

[0016] Signal masking refers to the process of embedding a message in a free-running chaotic time series such that the message is hidden, the chaotic dynamics are not altered enough to preclude synchronization, and the chaotic noise-like waveform is the dominant transmit feature. Signal masking consists of the message signal being added to a chaotic signal and recovered via subtraction of the synchronized receiver local chaotic signal from the received waveform. The masking signal requirements are satisfied when the message power has an embedding factor on the order of −30 dBc, which places the message power 30 dB below the chaotic power. Synchronization has been shown at higher information signal power levels, but with increased incidence of synchronization error bursts and degraded message recovery. One difficulty with this method is that noise power on the order of the message signal power level results in erroneous message recovery. The transmit SNR, which is determined by the power of the chaotic signal, must be very large in order to overcome the effects of modest channel noise which wastes power and broadcasts the signal presence (non-LPD). In addition, small trajectory perturbations, as little as the order of 10⁻¹⁶, can result in significant synchronization error bursts. Another difficulty for real implementations is that most synchronous chaotic systems can't accommodate channel losses, because the receiver synchronization degrades when the received signal has been attenuated from its original transmitted value. Finally, an unintended listener can recover the message without knowledge of the chaotic system by performing a mathematical reconstruction of the chaotic attractor and using it in conjunction with a set of recent observations to make local predictions of the immediate future for the chaotic system. Sufficiently accurate results are obtained to subtract out the chaotic process and uncover the hidden message, thereby compromising the security of the chaotic communication.

[0017] Chaotic attractors consist of an infinite number of unstable periodic orbits, all of which will be approached eventually from any initial condition. The system trajectory will be repelled away from any given orbit because of the stretching and folding of the chaotic dynamics, but it will eventually return to within a small neighborhood of the same orbit for the same reason. Orbit control exploits this type of behavior to transmit information by controlling the system trajectory such that it stays in the vicinity of one or more orbits for a prescribed length of time or follows a predetermined sequence of orbits. This method is very sensitive to noise and relatively small perturbations to the transmitted signal can cause the receiver trajectory to drift out of the required orbit.

[0018] Using the chaos shift key method of encoding digital data into chaotic transmitted sequences, for a binary-valued data stream, an aspect of the chaotic system is altered to create a transmit data stream. One method is to toggle a chaotic system parameter value, causing the data stream to switch between two chaotic attractors. Alternatively, the output can be inverted or not (i.e., multiplied by +1 or −1), causing the two attractors to be mirror images of each other. Both methods utilize a chaotic synchronization scheme as a precondition of message recovery. The message is demodulated either via the detection of synchronization with one of two attractors in the receiver (a hardware version of which is used with a Chua's circuit), or by the statistical method of correlation between the received sequence and a local chaotic signal over a large number of samples.

[0019] The orthogonal autocorrelation sequence method of encoding digital data into chaotic transmitted sequences develops a system with signal properties such that multiple user signals can be combined at the transmitter and then separated at each receiver by an autocorrelation calculation. The technique is presented as a general design method independent of the chaotic system used for the transmitter, but heavy use is made of a specific chaotic process in developing the concepts and results. Reasons for developing this scheme include the claims that non-repeating chaotic carrier signals provide greater transmission security than periodic carriers, require neither acquisition nor tracking logic for successful detection, and are more immune to co-channel interference than periodic carrier methods. Implementation investigations have demonstrated that an extremely large number of received values is required to accomplish demodulation, severely restricting the message data rates.

[0020] The present invention is a communications system based on a chaotic estimation engine. The invention provides for efficient communication with as little as two iterates per bit, a very general synchronization scheme, the integration of statistical methods with chaotic dynamics in the data recovery calculations, algorithms to enable a maximum aposteriori transmitted value estimation, secure communication by using completely overlapping ranges of transmit values for the conveyance of a logical zero or a logical one, and the development of a signal-to-noise ratio (SNR) estimator that produces good estimates of the true channel SNR and tracks changes in the channel noise power over time.

[0021] A transmit sequence is derived that sends one version of the chaotic strange attractor for a logical one and another version for a logical zero. Several new computational techniques allow for the estimation of both the instantaneous and average signal-to-noise ratio (SNR) values, as well as the maximum a posteriori (MAP) transmitted value for a given received chaotic iterate.

[0022] The receiver synchronizes to the transmitted chaotic sequence without requiring a chaotic equation separation into stable and unstable components or the design of an invertible chaotic process as was done in previous systems. The two processes of chaotic synchronization and message recovery required for synchronous communications are interleaved in the receiver estimation engine algorithms, obviating the need for an independent effort to achieve the synchronization of the transmitted chaotic signal with a locally generated version in the receiver. The estimation engine has the ability to make accurate decisions with as little as two chaotic iterations per data bit.

[0023] Three different estimates are generated for each received value and then combined into a single initial decision. This initial decision is mapped onto the closest portion of the geometric structures used to model the chaotic attractors via a minimum Euclidean metric. This mapping results in the final receiver decision of the transmitted value for the current received iterate. The concluding step of combining the decisions for all iterates in a bit yields the receiver bit decision, which is the transmitted message recovery.

[0024] A nonlinear chaotic receiver comprises a component for receiving a chaotic encoded digital signal transmission from a chaotic transmitter, synchronizing the chaotic receiver with the chaotic transmitter and recovering the contents of the encoded chaotic digital signal transmission using a chaotic strange attractor model and a chaotic probability density function model. Synchronization of the chaotic receiver with the chaotic transmitter and recovery of the contents of the encoded chaotic digital signal transmission may occur in the same calculations and result concurrently from the same calculations. The chaotic encoded digital signal transmission is a data sequence comprising a first through N number of iterates, wherein the first iterate represents a first value in the data sequence and the Nth iterate represents a last value in the data sequence. The chaotic strange attractor model comprises using one chaotic sequence to represent a logical zero data state and using a second chaotic sequence to represent a logical one data state. The first and second chaotic sequences may have attractors with completely overlapping regions of validity and the attractors may be mirror images of each other. The attractors are modeled as a set of geometrical functions having defined regions of validity. The chaotic strange attractor model comprises a strange attractor generated by combining Henon and mirrored Henon attractors, wherein the Henon and mirrored Henon attractors are generated by starting with one or more arbitrary points within an area of phase space that stretches and folds back onto itself, and inputting the points to a set of Henon equations, the result being the Henon attractor, and taking a mirror image of the Henon attractor to form the mirrored Henon attractor where the strange attractor is represented as a set of parabolas displayed on a Cartesian coordinate system and the parabolic regions of validity of the strange attractor are determined. The chaotic attractor model determines any existing fixed point on the strange attractor that repeats itself through multiple iterations of the chaotic transmission. The data sequence of the received chaotic encoded digital signal transmission is randomly selected from the group consisting of a first logical state for the Henon attractor and a second logical state for the mirrored Henon attractor. The chaotic probability density function models the probability of the first and second logical states of the Henon and mirrored Henon attractors as a random selection. The Henon strange attractor is generated by using image calculations on a Henon map, is represented in a Cartesian coordinate system as a crescent-like shape which occupies all four quadrants of the Cartesian coordinate system and is modeled as a set of four parabolas. The contents of the encoded chaotic digital signal transmission is determined by generating an initial decision for the contents of each iterate and generating a final decision for the contents of each iterate using a decision and weighting function. The final decision is generated using a synchronizer and final decision function and the final decision is used to recover the data sequence. The contents of the encoded chaotic digital signal transmission are determined by generating three estimates for each iterate received and calculating an initial decision for each iterate and mapping the initial decision onto the chaotic attractor to form a final decision for each estimate. The three estimates comprise a first estimate which is the value of the received iterate, a second estimate which is a minimum error probabilistic estimate and a third estimate which is a final decision of the previous iterate processed through Henon and mirrored Henon equations. The three estimates are combined to form the initial decision through a weighted average using probability calculations for the first, second and third estimates.

[0025] A means for determining a synchronization estimate to synchronize the chaotic receiver with a chaotic transmitter that generates the encoded chaotic digital signal transmission is provided. Signal estimates of the value of the 1 through N iterates received are generated.

[0026] The signal estimates comprise a received value which is equal to the actual value of the received iterate, a maximum a posteriori (MAP) estimate and a decision feedback estimate. The decision feedback estimate comprises generating a first decision feedback estimate for iterate n by passing iterate (n−1) through Henon equations and generating a second decision feedback estimate for iterate n by passing iterate (n−1) through mirrored Henon equations.

[0027] The decision and weighting function comprises performing a probabilistic combination of the three signal estimates by finding a probability for each iterate for each received value, MAP estimate and decision feedback estimate, generating an initial decision for each iterate, selecting the feedback estimate closest in value to the received value and the MAP estimate, generating weighting factors used in a weighted average computation and performing a weighted average, calculating a discount factor to discount a final decision for received values close to zero, determining the Cartesian coordinate representation of the initial decision and passing the Cartesian coordinates of the initial decision to a synchronizer and final decision generator and passing the zero proximity discount value to a system to bit map function. The weighting factors comprise generating a first weighting factor that quantifies channel effects by using a transmit PDF window used in determining the MAP estimate and multiplying the channel effects weighting factor by a second weighting factor that uses a PDF selected from the group consisting of the received PDF and transmit PDF, the result being the weighting factor for each iterate which is passed to the synchronization and final decision function. A discount weight is calculated to reduce the impact of initial decisions whose receive values are close to zero.

[0028] The selection of the feedback estimate comprises calculating one dimensional Euclidean distances from the received value and MAP estimate to the first and second feedback estimates and selecting the feedback estimate with the minimum Euclidean distance.

[0029] The synchronization and final decision function uses an initial decision from the decision and weighting function and maps Cartesian coordinates of the initial decision onto a Cartesian coordinate representation of the strange attractor. The mapping comprises inputting the initial decision into parabola models, the parabola models comprising equations approximating Henon and mirrored Henon attractors and calculating two dimensional Euclidean distance between the initial decision and every point on all parabolas and selecting a point on an attractor model corresponding to a minimum Euclidean distance as the final decision of the value of each iterate. The transmitted data sequence is recovered by combining the final decision of iterates 1 through N to determine the encoded digital signal transmission data sequence.

[0030] The nonlinear chaotic receiver comprises a receiver estimation engine for synchronizing the chaotic receiver with a chaotic transmitter and recovering the value of an encoded chaotic digital signal transmission. The receiver estimation engine comprises a signal-to-noise ratio (SNR) estimator, a maximum a posteriori (MAP) estimator, a feedback estimator wherein the chaotic receiver and the chaotic transmitter synchronization and the encoded digital signal transmission recovery occur concurrently while executing the same set of calculations within the receiver estimation engine. The receiver also includes a decision and weighting function within the receiver estimation engine comprising determining the probability of the estimate produced by the SNR estimator, the MAP estimator and the feedback estimator for each received iterate. The receiver calculates an initial decision for the iterate, a discount weight for a final decision for received values in close proximity to zero, and determines the final estimate of each iterate based on the initial decision from the decision and weighting function and then combines the final decision of iterates 1 through N to recover the encoded digital signal transmission data sequence.

[0031] The present invention comprises a method in a computer system and computer-readable instructions for controlling a computer system for receiving and recovering the contents of a chaotic encoded digital signal transmission by receiving a chaotic encoded digital signal transmission from a chaotic transmitter, synchronizing the chaotic receiver with the chaotic transmitter, and recovering the contents of the encoded chaotic digital signal transmission using a chaotic strange attractor model and a chaotic probability density function model.

[0032] These and other features, aspects and advantages of the present invention will become better understood with regard to the following description, appended claims and accompanying drawings where:

[0033] FIG. 1 shows a strange attractor of the Rossler system.

[0034] FIG. 2 shows a Poincare surface with intersecting trajectories.

[0035] FIG. 3 is a block diagram of a chaotic communication system.

[0036] FIG. 4 is a flowchart of the transmitter function.

[0037] FIG. 5 is a flowcart to model the chaotic attractor.

[0038] FIG. 6ais a plot of the result of generating a chaotic time series where an entire sector of the basin of attraction was chosen as the initial array condition.

[0039] FIG. 6bis a plot of the first iterate of the initial array of FIG. 6a.

[0040] FIG. 6cis a plot of the second iterate of the initial array of FIG. 6a.

[0041] FIG. 6dis a plot of the third iterate of the initial array of FIG. 6a.

[0042] FIG. 7 is a plot of a Henon attractor.

[0043] FIG. 8 is a plot of a mirrored Henon attractor.

[0044] FIG. 9 is a plot of the superposition of the attractors of FIGS. 7 and 8.

[0045] FIGS. 10athrough 10dare plots of parabolas modeling the Henon attractor. FIG. 10ais a plot of the innermost major section. FIG. 10bis a plot of the innermost minor section. FIG. 10cis a plot of the outermost minor section. FIG. 10dis a plot of the outermost major section.

[0046] FIG. 11ais a plot of the regions of validity for the Henon attractors of FIGS. 10athrough 10d.

[0047] FIG. 11bis a plot of the regions of validity of the parabolas of FIGS. 10athrough 10d.

[0048] FIGS. 12aand 12bare plots of examples of parabolas transposed through Henon equations.

[0049] FIG. 13 is a plot of parabolic transposition through the Henon map highlighting a fixed point

[0050] FIG. 14 is a plot of the transmit sequence empirical PDF for a data run of 4,000,000 iterates.

[0051] FIG. 15 is a plot of the transmit PDF of FIG. 14 overlaid on the PDF model.

[0052] FIG. 16 is a block diagram of the receiver and the system to bit map function.

[0053] FIG. 17 is a flowchart of the receiver estimation engine process.

[0054] FIG. 18 is a flowchart of the SNR maximum probability estimator.

[0055] FIG. 19aare plots of the PDF traces for an SNR of 6 dB.

[0056] FIG. 19bare plots of the PDF traces for an SNR of 15 dB.

[0057] FIG. 20aare plots of the PDF traces for an SNR of 30 dB.

[0058] FIG. 20bare plots of the PDF traces for an SNR of 60 dB.

[0059] FIG. 21 is a plot of a two dimensional probability function at about 30° rotation and 120° elevation.

[0060] FIG. 22 is a plot of a two dimensional probability function at about 5° rotation and 10° elevation.

[0061] FIG. 23 shows a plot of PDF dependency on SNR for various received iterate values.

[0062] FIG. 24 is a plot of the SNR estimator performance for a 64 iterate window.

[0063] FIG. 25 shows a flowchart of the general MAP transmit value determination where the transmit sequence PDF is determined and stored by empirical observation or using model calculations at receiver initialization.

[0064] FIG. 26 is a flowchart of the MAP transmit value determination with Gaussian noise.

[0065] FIG. 27 is a flowchart of the simplified MAP value determination.

[0066] FIG. 28 is a graphical representation of the transmit PDF window construction.

[0067] FIG. 29 is a graphical representation of the MAP decision windowing and maximum value determination.

[0068] FIG. 30 is a graphical representation of the probability of reversed noise PDF centered on the received iterate and evaluated at zero.

[0069] FIG. 31 is a flowchart of the synchronizer/final decision function.

[0070] FIG. 32 is a plot of the results of the receiver estimation engine model for SNR values between 0 and 26 dB and for a number of iterates per bit of 2, 3, 4, 8, 16, 32 and 64.

[0071] FIG. 1 is a depiction of the strange attractor of the Rossler system, showing the chaotic characteristics of stretching and folding. As a background, chaotic dynamics are frequently viewed in phase space, which is a region of space where the current state of the system perpetually exists. Stretching 10, has causes such as scaling factors greater than unity, power terms such as squared, cubed, etc., or other mathematical processes. The stretching expands the phase space and causes nearby points 11 and 12 in the phase space to move apart as time progresses. Folding 13 is a result of nonlinearities such as power terms or appropriate slope breakpoints. Folding causes the stretched phase space to fold back into its original confinements, resulting in the capture of the chaotic dynamics in a finite area of space 14. Phase space and the structure contained within can be illustrated with a chaotic map (i.e., a set of equations with chaotic behavior) in three-dimensional space using the Rossler system. The chaotic equations operate on coordinates in space, and each iteration of the equations signifies the passage of the next time increment.

[0072] The Rossler equations in differential form are $\begin{array}{cc}\frac{dx}{dt}=-\left(y\left(t\right)+z\left(t\right)\right)\frac{dy}{dt}=x\left(t\right)+\mathrm{ay}\left(t\right)\frac{dz}{dt}=b+x\left(t\right)z\left(t\right)-\mathrm{cz}\left(t\right)& \left(1\right)\end{array}$

[0073] where a=0.2, b=0.2, and c=5.7 in this case, and t denotes time. In discrete time form, these equations become

[0074] x[k]=x[k−1]−(y[k−1]+z[k−1])

[0075] y[k]=y[k−1]+x[k−1]+ay[k−1]

[0076] z[k]=z[k−1]+b+x[k−1]z[k−1]−cz[k−1]  (2)

[0077] where k is the discrete time index. For any time k, the coordinate value (x[k],y[k],z[k]) is plotted in space. The collection of points for a large number of iterations in the region of space where the Rossler chaotic system perpetually exists is shown in FIG. 1. The structure that forms from the continued iteration of the chaotic map in this region of perpetual chaotic existence is called the strange attractor of the chaotic system. The strange attractor can be visualized in different ways. One is simply the collection of points that results from repeated iterations of the chaotic system, neither requiring knowledge of the temporal sequence of the points nor using connecting lines between the data points. It is the purest form of the attractor and applies to all attractors, but cannot illustrate the dynamics associated with the chaotic process time evolution. Another attractor visualization method uses a collection of trajectories or orbits that the chaotic system follows. A trajectory or orbit is a sequence of points produced by the dynamics of the chaotic system. Its relationship to the strange attractor is best illustrated by plotting the points in phase space in their temporal order of occurrence, with lines connecting them. The sequence of phase space location can then be traced on the attractor, and such processes as stretching and folding might be illustrated. Some attractors are better than others for this method of visualization. FIG. 1 is an illustration of this trajectory plotting method of attractor construction.

[0078] Under iteration, these systems induce all points in phase space within a basin of attraction to collapse to a chaotic attractor and all points outside the basin of attraction to fly off to infinity. The attractor is a fractal, which is a geometrical figure with non-integer dimension. It is comprised of an infinite number of unstable orbits.

[0079] A useful viewpoint of a chaotic process is found by creating a cross-section of its orbit at some position in phase space as shown in FIG. 2, which is a general illustration of a Poincare surface with intersecting trajectories. FIG. 2 shows this cross-section 20, called a Poincare section, and the pattern of intersections for the trajectory crossing the surface of section 21 can have some useful applications in controlling the orbits of a chaotic system either to mitigate chaotic effects in applications where chaos is undesirable or to impart a property to the chaotic system such as information content in a communications system. Every time that a chaotic trajectory 21 crosses the surface of section 20, its position 23 on the section can be determined and the orbit modified if desired. Orbit modification may minimize the drift inherent in one of the infinite number of unstable orbits contained in the chaotic attractor so a predetermined orbit is perpetually maintained. Another possibility is to alter the orbit so it follows a particular sequence of states in the trajectory. A third approach may be to eliminate states from the orbit for the purpose of encoding information into a transmitted chaotic sequence.

[0080] The top-level chaotic communications system block diagram is depicted in FIG. 3, showing a transmitter 25, channel 26, which may be a lossless non-fading AWGN, and receiver 27. In the transmitter, the chaotic system is specified as a Henon map 28, and the alternate version of the same strange attractor labeled the mirrored Henon map 29. Typical data sequences have a stream of logical 0 and 1 values, each occurring with a probability of ½. In these simulations, the data bits are therefore generated as a random bit stream 30, and the logical values of 0 and 1 toggle a switch 31 between the Henon and mirrored Henon chaotic systems.

[0081] Although not explicitly shown in FIG. 3, the output of the switch 31 is passed through a digital-to-analog converter (D/A), followed by a power amplifier and possibly frequency translation in order to transmit the encoded chaotic sequence. The practical issues of bit alignment, sampling synchronization, D/A and A/D quantization effects, and channel phenomena such as losses, fading, Doppler, and non-Gaussian noise are addressed in the hardware implementation of this chaotic communications technique.

[0082] The receiver 27 in FIG. 3 shows the estimation engine 32, followed by the data bit calculation block 33. Transmitted value and transmitted chaotic system decisions are made on every received value. Data bits containing a certain number of chaotic iterates are recovered by combining the appropriate individual iterate decisions to reach a bit decision to decode the chaotic transmission.

[0083] The transmitter 25 of FIG. 3 is based on the Henon chaotic system. FIG. 4 shows a flowchart of the transmitter function 40. Before transmitting the encoded chaotic signal 45, the transmit function must operate free-running chaotic oscillators 41 and 42, generate the transmission data bit stream 43, encode the data as chaotic iterates 44, convert the iterates from digital to analog and amplify and frequency translate the signal if necessary 45. In designing the receiver estimation engine, the strange attractor and the probability density function must be empirically determined in any implementation of the receiver estimation algorithms. In addition, either or both may be modeled for certain algorithm designs.

[0084] The modeling of the chaotic attractor may proceed as shown in FIG. 5. The first step in modeling the chaotic attractor 50 is to generate the Henon attractor 51 and the mirrored Henon attractor 52. The next step is to create models of the strange attractor 53. After the strange attractor is modeled 53, the next steps are to determine the parabolic regions of validity 54, the chaotic dynamics on the parabolic regions 55 and to determine the Henon fixed point 56. The processing of these steps is detailed below.

[0085] A standard attractor generation technique involves creating a chaotic time series by starting at a single initial point and iterating through the chaotic equations until sufficient information has been obtained to visualize an attractor. This process is shown in FIGS. 6athrough d. Instead of starting with a single point and generating a plot from the time series, an entire section of the basin of attraction was chosen as the initial array condition (FIG. 6a). This region quickly changed its morphological appearance into a form resembling the Henon attractor. FIGS. 6b, 6c, and 6dshow the first, second, and third iterates of the initial array respectively. The choice of a rectangular subset of the basin of attraction for the process as shown in FIGS. 6athrough 6dwas one of convenience. Any other area of arbitrary shape wholly contained within the basin of attraction would have generated the same result, with more or fewer iterations necessary, depending on its size.

[0086] FIG. 7 depicts the resulting Henon attractor. The equations for the Henon chaotic process are that produced this attractor are:

[0087] x_h[k+1]=−ax_h[k]²+y_h[k]+1  (3)

[0088] y_h[k+1]=bx_h[k]k=1,K,K  (4)

[0089] where a=1.4, b=0.3, k=iterate index, and there are K total iterates in a time series.

[0090] This process for generating the Henon attractor works because there is a trapping region for chaotic systems, also called a basin of attraction, within which no orbit can escape. In general space, the stretching and folding of the chaotic dynamics will cause many orbits to escape to infinity, but there exists a region of phase space which stretches and folds back onto itself ad infinitum. This region of phase space is called the basin of attraction. This entire trapping region maps back onto itself with a morphologically altered appearance, and multiple iterations of the chaotic transform result in the strange attractor. Letting R represent the basin of attraction and H represent the Henon transform,

[0091] H(R)⊂R,  (5)

[0092] so the Henon transformation of R lies entirely within R. Repeated applications of this region produce subsets of the region, and so no orbits can escape.

[0093] FIG. 7 is a plot of a Henon attractor and is the superposition of iterates 5-10,000 of this initial set. The first few iterates of the Henon transform were eliminated because they contained excess area tending to obscure the attractor detail. Iterates in the vicinity of 10,000 had only a few points left to plot because the finite precision of the internal representation of numbers inside computers caused the transforms of distinct numbers to be indistinguishable to the computer. The range of x-values for the Henon system was approximately symmetric, covering roughly (−1.28, 1.28). The y-values have the same approximate symmetry over a smaller range because the they are scaled versions of the x-values, as can be seen in equation (4). This symmetry allows the creation of a second version of the same attractor, mirrored through the x- and y-axes. This mirrored attractor has an identical range of x-values, yielding the desired transmit characteristic of completely overlapping transmit values for both logical data states. It was observed that the structure of the main attractor features appeared roughly parabolic. Simple geometric equations can, therefore, be used to model the major attractor sections, yielding a means by which to map a probability-based receiver decision onto the chaotic dynamics represented by the attractor.

[0094] Turning back to FIG. 5, the next step is to empirically determine the mirrored Henon attractor 52 using an identical process as that depicted in FIG. 6 for the Henon chaotic process. The output of that step is a mirrored Henon attractor, shown in FIG. 8, which is a second version of the attractor shown in FIG. 7 created for the second logical state. It is a mirror image through the x- and y-axes which yields intersecting crescents with completely overlapping x-value ranges. Letting

[0095] x_m=−x_hand y_m=−y_h,  (6)

[0096] x_m[k+1]=1.4x_m[k]²+y_m[k]−1  (7)

[0097] y_m[k+1]=0.3x_m[k].  (8)

[0098] The attractor formed by these equations is shown in FIG. 8, along with the corresponding basin of attraction. FIG. 9 shows the superposition of the two versions of the strange attractor (FIGS. 7 and 8) and highlights the identical ranges of values for both the Henon and mirrored Henon systems. The complete overlap of value ranges is indicative of the fact that there is not a one-to-one correspondence between a transmit value or a range of transmit values and digital message value, as exists for most other communications systems. This feature yields security for the communications scheme because the observation of a received value or range of values does not reveal in a direct manner which logical data state has been sent. The sequence of received values in a data bit must be processed and the attractor from which the sequence was generated found in order to determine the logical data state of the bit being received. There are areas where the two attractors intersect each other 70 and 71, and it will be difficult to distinguish between the Henon or mirrored Henon chaotic dynamics in these regions. In addition, there are precise points of intersection for which the correct decision cannot be guaranteed even in the absence of noise.

[0099] Turning back to FIG. 5, the next step is to model the Henon strange attractor 53. The Henon strange attractor 53 may be modeled with four parabolas to achieve a fairly accurate representation of the structure, and the plots are shown in FIGS. 10athrough 10d. Image tools were used again to plot the attractor and parabolas on the same image window to evaluate the match of a given parabola to the appropriate attractor section. The Henon attractor was broken into a sections according to a large (major) crescent pair, with a smaller (minor) crescent pair contained inside the major section. These four primary features were deemed to model the attractor well, and so four parabolas were found. The parabola equations and valid y-values are listed below in Table 1 for each of the plots of FIGS. 10athrough d.

[0100] To quantify the closeness of the parabolic model to the empirical attractor, a time sequence of Heron (x,y) values was generated and the y-values used as input to the four parabola equations. This created four corresponding x-values for each Henon x-value, and the closest parabola was chosen for each Henon point. The resulting set of errors was obtained by subtraction of only the x-values for each parabola (because the Henon y-values were used in parabola equations), and the mean-square-error (mse) was calculated for each parabola fit. Increasingly longer time sequences were used until the results stabilized. A 500,000 point time sequence generated the results shown in Table 2, which illustrates the excellent quality of fit between the empirically determined attractors and their parabolic section models.

[0101] Turning back to FIG. 5, the next step in modeling the chaotic attractor 50 is to determine the parabolic regions of validity 54. The regions of validity for the four parabolas may be determined by the limits of the major (large pair of crescents) and minor (small pair of crescents) sections on the attractor. These regions are shown in FIGS. 11aand 11b, where FIG. 11ashows the Henon attractor and FIG. 11bshows the four modeling parabolas.

[0102] Referring to FIG. 11a, the vertical dash marks 73 delineate the limits of the major and minor sections on the attractor. These limit marks 73 are transposed directly onto the parabolas in FIG. 11b, identifying almost all termination points on the parabolic regions of validity. There is, however, an arrow 74 in the vicinity of point (0.6,−0.22). This arrow 74 marks the limit of the parabola 3 region of validity, because this area of the attractor is about equidistant between parabolas 3 and 4, and the model transitions to parabola 4 at this point.

[0103] Turning back to FIG. 5, the next step in modeling the chaotic attractor 50 is to determine the chaotic dynamics associated with the various parabolic regions 55. Each of the model parabolas was broken into sections according to quadrant of the Cartesian coordinate system and traced through an iteration of the Henon system equations to determine which parabolic sections they became after a chaotic iteration. Examples are shown in FIGS. 12aand 12b. Each final receiver decision of the current transmit value is mapped by a minimum Euclidean distance metric onto the closest attractor section via the parabolic model equations. Assuming an accurate decision, knowledge of the current parabolic section and its transposition through the system equations can be used to identify the range of valid received values on the next iterate. This information is embedded in other areas of the receiver calculations, and so these tracings are no longer explicitly used. However, while determining the chaotic dynamics associated with the various parabolic regions, a fixed point on the Henon attractor was identified.

[0104] Turning back to FIG. 5, the last step in modeling the chaotic attractor 50 is to determine the Henon fixed point 56. Tracing the parabolic sections through the Henon map uncovered a repelling fixed point 75 in the Henon attractor which is shown in FIG. 13. A fixed point on a chaotic attractor is one that repeats itself through chaotic iterations ad infinitum, i.e. it maps exactly onto itself. A stable or attractive fixed point is one toward which a neighborhood around the fixed moves under iteration of the chaotic equations, so chaotic iterates on the neighborhood decrease the distance from the fixed point. An unstable or repelling fixed point 75 is one for which the distance between itself and any point in a surrounding neighborhood increases through iteration of the chaotic system. Within the resolution of MATLAB, the Henon fixed point is (0.63135447708950, 0.18940634312685). This point has been found to the 14^thdecimal place. The first iterate changes the x-value by one digit in the 14^thdecimal place, with the y-value unchanged. After 45 iterations of the Henon equations, the result is still within 2% of the starting values.

[0105] The chaotic attractor has been empirically determined and modeled as shown in FIG. 5. The next step in the determination of a priori information for proper receiver design is to find, and possibly model the chaotic transmit sequence probability density function. The transmit sequence is chosen between the Henon and mirrored Henon free-running maps by the logical value of the data, with a logical 0 choosing the mirrored Henon value and a logical 1 choosing the Henon value for this investigation. A typical digital data stream contains about equal occurrences of zeros and ones, so the transmit data sequence was modeled as a random sequence of zeros and ones with equal probability of occurrence. The PDFs for the Henon, mirrored Henon, and transmit sequences were found empirically using multimillion iterate data runs. All three PDFs were modeled, but the transmit PDF is the only used in the receiver estimation engine for this design.

[0106] Starting from any given initial point within the appropriate basin of attraction, the first 100 or so iterates were ignored as transient while the chaotic trajectory moved onto the attractor. A time series was then recorded for 2,000,000 iterates for the Henon and mirrored Henon systems, and 4,000,000 iterates for the transmit sequence.

[0107] It was desirable to define the range of values according to the anticipated values for the received data sequence with noise contamination. This range of values was based on a worst-case signal-to-noise ratio (SNR) of 0 dB, at which point the noise power and signal power are equal. It is generally difficult to work with a higher noise power than signal power, and most techniques approach a 0.5 BER in these cases. It was decided to use 4096 data bins to perform the various PDF calculations and the subsequent receiver processing on all of the PDFs so found.

[0108] As discussed below in more detail concerning the receiver processing, the received data PDF is the convolution of the transmit PDF and the PDF of the additive noise in the channel. In the general case of a received PDF without closed form representation, the FFT-based fast convolution technique can be used to find the received PDF for any SNR value and any combination of transmit and noise PDF waveforms. It is therefore important to define the bin range with a sufficient region of insignificant values to make the aliasing contamination of the desired linear convolution result caused by the inherent circular convolutional characteristic of the FFT-based method negligible. In addition, the greatest computational efficiency of the FFT is realized when its input data length is a power of 2. Zero-padding can be used to address the circular convolution concern. However, the increased lengths of the FFT calculations load down the computer system resources if 4096 bins are spread over a smaller range of received values such that zero padding is required, and the calculation fidelity decreases if fewer bins are used over the range of significant values. A range of [−5.12, 5.12] with 4096 bins was used, yielding a bin width of $1/400,$

[0109] or 0.0025, and a bin density of $400\mathrm{bins}/\mathrm{UnitRange}.$

[0110] To empirically find the many PDFs in this investigation, the MATLAB histogram function was used for each PDF calculation to get the number of occurrences in the time series that fell into the range of numbers for each data bin.

[0111] Letting b=bin index

[0112] B=number of data bins

[0113] K=number of iterates in the time series

[0114] _bn=number of hits in data bin b

[0115] _br=relative number of occurrences in data bin b

[0116] w=bin width, which is constant for all bins

[0117] p=continuous probability density function

[0118] _bp=PDF value of data bin b (quantized probability density function)

[0119] P=probability $\begin{array}{cc}{\hspace{0.17em}}_{b}r={\hspace{0.17em}}_{b}n/K& \left(17\right)\end{array}$

[0120] where $\begin{array}{cc}\sum _{b=1}^B{\hspace{0.17em}}_{b}n=K& \left(18\right)\end{array}$

[0121] and _bn is found easily by using MATLAB's histogram function.

[0122] The fundamental definition of PDF requires that $\begin{array}{cc}{P}_{\infty }={\int }_{-\infty }^{\infty }\left(p\left(x\right)\right)dx=1.& \left(19\right)\end{array}$

[0123] For the present formulation, this translates into the requirement $\begin{array}{cc}\sum _{b=1}^B{\hspace{0.17em}}_b\mathrm{pw}=1.& \left(20\right)\end{array}$

[0124] The PDF can be found from the relative number of occurrences as follows: $\begin{array}{cc}\sum _{b=1}^B{\hspace{0.17em}}_{b}r=\sum _{b=1}^B{\hspace{0.17em}}_{b}n/K=K/K=1& \left(21\right)\\ \sum _{b=1}^B{\hspace{0.17em}}_b\mathrm{pw}=w\sum _{b=1}^B{\hspace{0.17em}}_{b}p=\sum _{b=1}^B{\hspace{0.17em}}_{b}r& \left(22\right)\\ {\hspace{0.17em}}_{b}p={\hspace{0.17em}}_{b}r/w.& \left(23\right)\end{array}$

[0125] Because the logical 0 and 1 values of a typical data sequence each have a probability of occurrence of ½, the transmit sequence of the transmit probability density function is modeled as an arbitrary choice between the Henon and mirrored Henon time series for each data bit duration. Data runs of 4,000,000 iterates were taken to empirically construct a smooth transmit PDF as shown in FIG. 14. This transmit PDF can be found from the Henon and mirrored Henon PDFs via multiplying each by its corresponding probability of occurrence (½ for both in this case) and adding them together. Actually, all PDFs can be generated from a single Henon or mirrored Henon PDF. The transmit PDF was modeled as a DC level plus a summation of Gaussian functions. There is a choice as to window only the DC level and let the Gaussian functions remain infinite, or to set up the DC level and Gaussian functions first and window the entire PDF model. In either case, the window must be defined from observations of the empirical PDF trace coupled with the bin range and bin width choices. The maximum positive bin center value for which there is a nonzero transmit PDF value is 1.2841, and the bin limits are [1.2828, 1.2854]. Since the transmit PDF is symmetrical, the corresponding negative bin limits are [−1.2828, −1.2854]. The window is, therefore, defined as

[0126] W_dc(x)=U(x−(−H_max))−U(x−H_max)=U(x−(−1.2854))−U(x−1.2854)  (24)

[0127] where $U\left(x_0\right)=\left\{\begin{array}{c}0\Rightarrow x<x_0\\ 0\Rightarrow x\ge x_0\end{array}\right\}=\mathrm{the}\mathrm{unit}\mathrm{step}\mathrm{function}$

[0128] and H_max=1.2854 is the maximum Henon magnitude for the specified bin structure

[0129] The differences between the two windowing methods are very small in terms of establishing the PDF trace model. It turns out that windowing both the DC level and the Gaussian functions requires the DC level (0.2160) to be greater by 0.0002 than the method of windowing only the DC value (0.2158). The small degree of difference is because the Gaussian functions used to model the transmit PDF have very little area outside the valid DC region. The choice between windowing methods cannot be made on considerations stemming from the construction of a PDF model, but SNR estimation calculations involving the transmit PDF in the receiver influenced the decision to window only the DC component.

[0130] First the DC level and four Gaussian functions structure the general trend with MSE (mean square error)=0.0057. Some coarse details are added in with 18 more Gaussian functions (22 total) with MSE=0.0031, and some detailed fine structure consisting of 16 Gaussian functions (38 total) and MSE=0.0014. The final attractor approximation has the final details added in using 6 Gaussian functions (44 total) that fill in some gaps and have a final MSE=0.0012.

[0131] The final tally is a windowed DC level with 44 weighted Gaussian functions to construct the transmit PDF model. Because of the symmetry of the structure, only 22 Gaussian functions (mean and variance) had to be defined, along with the corresponding 22 weighting factors.

[0132] The empirical attractor of the transmit PDF is overlaid on the model as shown in FIG. 15. The close match with the distinguishing features of the true PDF can be readily seen. A metric against which to measure the match between the empirical PDF and the approximation is the mean squared error (MSE), whose value is 0.0012. The approximation total probability, which is the area under the PDF approximation, is 0.9999. This compares favorably with the theoretical area under a PDF curve, which equals one.

[0133] The equation for the transmit PDF approximation model with windowing only for the DC level is $\begin{array}{cc}{\tilde{p}}_t\left(x\right)={\mathrm{cW}}_{dc}\left(x\right)+\sum _{l=1}^L\frac{b_l}{\sqrt{2\pi }{\sigma }_l}{e}^{-{\left(\frac{x-m_l}{\sqrt{2}{\sigma }_l}\right)}^2}& \left(25\right)\end{array}$

[0134] where c=DC value

[0135] l=Gaussian functions' index

[0136] L=total number of Gaussian functions

[0137] b_l=weights of the Gaussian functions

[0138] m_l=mean values of the Gaussian functions

[0139] σ_l²=variances of the Gaussian functions

[0140] and the values for the approximation depicted here are listed below, with a one-to-one correspondence between the values occupying identical positions in the vectors b,m, and σ

[0141] c=0.2158

[0142] L=44

[0143] The power in the encoded chaotic transmit sequence needs to be determined. This can be done using the power found empirically in the Henon and mirrored Henon time signals. It should be noted that the Henon and mirrored Henon PDFs were skewed positive and negative, respectively, indicating correspondingly positive and negative mean values for the systems. The mean and variance for the systems have been estimated from million-iterate data runs to be

[0144] As shown in FIG. 15, the transmit PDF is seen to be symmetrical about zero, however, because of the random selection between the Henon and mirrored Henon time sequences according to the data bits. Since the Henon and mirrored Henon values ranges have complete overlap, the range of transmit values is the same as the range in either individual system. In terms of the original systems, $\begin{array}{cc}{\stackrel{\_}{x}}_h=\frac{1}{K}\sum _{k=1}^K{x}_h\left[k\right]& \left(26\right)\\ {\sigma }_h^2=\frac{1}{K}\sum _{k=1}^K{\left(x_h\left[k\right]-{\stackrel{\_}{x}}_h\right)}^2=\frac{1}{K}\sum _{k=1}^K{\stackrel{\_}{x}}_h^2\left[k\right]-{\stackrel{\_}{x}}_m^2=\stackrel{\_}{x_h^2}-{\stackrel{\_}{x}}_h^2& \left(27\right)\\ {\stackrel{\_}{x}}_m=\frac{1}{K}\sum _{k=1}^K{x}_m\left[k\right]& \left(28\right)\\ {\sigma }_m^2=\frac{1}{K}\sum _{k=1}^K{\left(x_m\left[k\right]-{\stackrel{\_}{x}}_m\right)}^2=\frac{1}{K}\sum _{k=1}^K{x}_m^2\left[k\right]-{\stackrel{\_}{x}}_m^2=\stackrel{\_}{x_m^2}-{\stackrel{\_}{x}}_m^2.& \left(29\right)\end{array}$

[0145] Turning back to FIG. 4, the encoded chaotic transmit sequence is now generated 44, modeled as a random selection between the Henon and mirrored Henon systems 43. For the transmit sequence,

[0146] P_h(x_h)=P_r(x_his sent) P_m(x_m)=P_r(x_mis sent)  (30)

[0147] P_h(x_h)=P_m(x_m)=0.5  (31)

[0148] x_t=(x_h)(P_h(x_h))+(x_m)(P_m(x_m))  (32)

[0149] $\begin{array}{cc}{\stackrel{\_}{x}}_t=\frac{1}{K}\sum _{k=1}^K{x}_t\left[k\right]=P_h\left(x_h\right)\left(\frac{1}{K}\sum _{k=1}^K{x}_h\left[k\right]\right)+P_m\left(x_m\right)\left(\frac{1}{K}\sum _{k=1}^K{x}_m\left[k\right]\right)& \left(33\right)\\ {\stackrel{\_}{x}}_t=\frac{{\stackrel{\_}{x}}_h+{\stackrel{\_}{x}}_m}{2}=\frac{0.2558-0.2558}{2}=0& \left(34\right)\\ {\sigma }_t^2=\frac{1}{K}\sum _{k=1}^K{\left(x_t\left[k\right]-0\right)}^2=P_h\left(x_h\right)\left(\frac{1}{K}\sum _{k=1}^K{x}_h^2\left[k\right]\right)+P_m\left(x_m\right)\left(\frac{1}{K}\sum _{k=1}^K{x}_m^2\left[k\right]\right).& \left(35\right)\end{array}$

[0150] Substituting from equations (27) and (29) into (35) yields the transmit power equation $\begin{array}{cc}{\sigma }_t^2=P_h\left(x_h\right)\left({\sigma }_h^2+{\stackrel{\_}{x}}_h^2\right)+P_m\left(x_m\right)\left({\sigma }_m^2+{\stackrel{\_}{x}}_m^2\right)=\frac{\left({\sigma }_h^2+{\stackrel{\_}{x}}_h^2\right)+\left({\sigma }_m^2+{\stackrel{\_}{x}}_m^2\right)}{2}.& \left(36\right)\end{array}$

[0151] and entering the empirically determined Henon and mirrored Henon values yields $\begin{array}{cc}{\sigma }_t^2=\frac{\left(0.5210+{0.2558}^2\right)+\left(0.5210+{\left(-0.2558\right)}^2\right)}{2}=0.5864.& \left(37\right)\end{array}$

[0152] {overscore (x)}_t=0

[0153] σ_t²=0.5864

[0154] These values have been proven by statistical analysis of empirical data runs. The transmit power is used to set up the proper lossless non-fading AWGN channel noise power level via an SNR specified as an input to the simulation for the bit error rate (BER) data runs.

[0155] The general block diagram for the receiver is shown in FIG. 16. When a chaotic transmission is received 100, three estimates are generated for each received iterate 100, an initial decision is calculated from them, and a final decision is formed by mapping the initial decision onto the chaotic attractor. The mapping operation also yields the synchronization estimate. It provides chaotic synchronization without requiring the special characteristics or techniques utilized by other methods, such as stable/unstable chaotic subsystem separation or inverse system methods. At the appropriate times, the individual iterate decisions are combined into a bit decision. A BER is then calculated to evaluate the performance of this method.

[0156] The first estimate is the received value 100 itself, called the received iterate. In the absence of noise, it is an accurate indication of the transmitted value. In the presence of additive white Gaussian noise, the received value is a valid estimate of the transmitted iterate because the most likely noise value is 0. The maximum a posteriori (MAP) transmit value 101 function calculates the second estimate for every received iterate and yields the minimal error probabilistic estimate. The Synchronization/Final Decision Generator 102 calculates the third estimate using the receiver final decision of the previous iterate, which is delayed by one iterate 107 and then processed through the Henon and mirrored Henon transform 103 equations. The third estimate provides information on the next possible received values for both chaotic systems given an accurate current decision, which aids the decision process in the presence of noise.

[0157] The Decision and Weighting function 104 combines the three estimates into an initial decision through a weighted average whose weights derive from probability calculations on each estimate. It is given the received iterate estimate, the MAP estimate, and two possible feedback estimates—one assuming a correct previous decision and the other assuming an erroneous previous decision. It chooses the feedback estimate closest to the received value and MAP estimates, determines the necessary probability values, and performs the weighted average calculation.

[0158] The Decision and Weighting function 104 also generates a discount weight to be used by the System to Bit Map function 105 in combining individual iterate decisions into a bit decision that serves to de-emphasize iterate decisions whose received value is in close proximity to zero. The Henon and mirrored Henon transmit possibilities are the negatives of each other. Transmitted values in the close proximity of zero can appear to have come from the incorrect attractor with very small amounts of noise. The feedback estimate chosen as closest to the MAP and received value estimates can easily be in error, which may then create an erroneous final iterate decision. The quantification of the term “close proximity” is based on the average SNR estimate 106, discounting more severely those received values close enough to zero to have been significantly affected by the noise currently present in the channel, while emphasizing the larger magnitude received values.

[0159] The Synchronization and Final Decision Generator 102 is given the initial probabilistic decision and the corresponding y-value from the Henon or mirrored Henon transform. It uses parabola equations that model the Henon and mirrored Henon attractors to generate parabolic sections in the vicinity of the input point and chooses the parabola and corresponding system closest to the input point via a Euclidean distance metric. This process adjusts both the x- and y-values of the receiver initial decision to provide accurate tracking of the chaotic dynamics in the final decision.

[0160] FIG. 17 is a flowchart of the receiver estimation engine process 110. As stated above, the three estimates generated for each received iterate are the received value itself, the MAP estimate, and the decision feedback estimate. The received value is used directly because it equals the transmitted value in the absence of noise and, because the AWGN has zero mean, the expected value of the received value is the transmitted value. The other two estimates must be generated by the receiver estimation engine as shown in FIG. 17. Upon initialization 111, the signal-to-noise (SNR) estimator parameters, chaotic feedback initial conditions, strange attractor model initial conditions (from the a priori receiver design investigation) and transmit PDF model calculations (from the a priori receiver design investigation) are input. The next step is to receive the transmitted chaotic signal and to process the first iterate received 112.

[0161] The next step is to estimate the average SNR 113. The MAP estimator requires knowledge of the noise power (or SNR). The SNR or noise power is also used to find receiver probabilities needed in the weighted averaging of the three estimates into an initial decision and, again, as the means for discounting iterate decisions that are too close to zero in the bit decision process. The receiver requires a method for estimating the noise power contained in a chaotic signal contaminated with random noise. The average SNR estimator 113 consists of an SNR maximum probability estimator 114 followed by a running weighted average calculation 115.

[0162] Turning now to FIG. 18, a flowchart is shown of the SNR maximum probability estimator 150 (114 in FIG. 17). The signal transmission is a data sequence comprising N number of iterates where the first iterate represents a first value in the data sequence and the Nth iterate represents a last value in the data sequence. For a received signal iterate value 151, a lookup table can be used to find the maximum probability SNR for a received iterate 152. This involves determining the instantaneous SNR, which can be defined as the instantaneous SNR for which a received PDF trace has the greatest probability of occurrence out of the PDF traces for all SNR values. The lookup table contains received signal PDFs that have been stored in memory 153. The received signal iterate value is passed into the lookup table and the corresponding PDF value is returned 154 for all SNR values. The SNR having the largest PDF value associated with the received value is chosen 155. The individual iterate decisions are then smoothed through a running average calculation 160 and processing continues in FIG. 17.

[0163] Alternatively, the instantaneous SNR may be calculated rather than using a lookup table 156. An algorithm is used and is based on the closed-form equations realized for the received PDF. The model chosen for the transmit PDF is based on Gaussian functions with an AWGN channel model 157. The received PDF is the convolution of the transmit PDF with the channel noise PDF 158, and the convolution of a Gaussian function with another Gaussian function is a Gaussian function. This closed-form representation enabled the development of the computational algorithm for the determination of the instantaneous maximum probability SNR.

[0164] The transmit PDF model that windows only the DC component and not the Gaussian functions is used. Although both transmit PDF models result in closed-form expressions for the received PDF, the version that windowed all terms had significantly more complicated expressions involving the evaluation of Q-functions in all terms plus exponential functions in most terms. The version that windowed only the DC component had two Q-functions for the windowed DC term, with the transmit PDF Gaussian terms remaining Gaussian in the received PDF.

[0165] The calculation of Q-functions is computationally more expensive than the evaluation of Gaussian functions, even when a Q-function approximation formula is used. This additional complexity is easily seen by comparing the general Gaussian equation in equation (38) with the integrated Gaussian structure of the Q-function definition equation (44) and the extra terms preceding the Gaussian function in the Q-function approximation equation (45).

[0166] A few terms will be collected from previous sections to continue the discussion. The equation for the transmit PDF approximation model (FIG. 18, 157) is $\begin{array}{cc}{\tilde{p}}_t\left(x\right)={\mathrm{cW}}_{dc}\left(x\right)+\sum _{l=1}^L\frac{b_l}{\sqrt{2\pi }{\sigma }_l}{e}^{-{\left(\frac{x-m_l}{\sqrt{2}{\sigma }_l}\right)}^2}& \left(25\right)\end{array}$

[0167] with

[0168] W_dc(x)=U(x−(−H_max))−U(x−H_max)=U(x−(−1.2854))−U(x−1.2854)  (24)

[0169] where $U\left(x-x_0\right)=\left\{\begin{array}{c}0\Rightarrow x<x_0\\ 1\Rightarrow x\ge x_0\end{array}\right\}=\mathrm{the}\mathrm{unit}\mathrm{step}\mathrm{function}$

[0170] and H_max=1.2854 is the maximum Henon magnitude for the specified bin structure. The channel noise PDF is modeled as $\begin{array}{cc}{p}_n\left(x\right)=\frac{1}{\sqrt{2\pi }{\sigma }_n}{e}^{-{\left(\frac{x}{\sqrt{2}{\sigma }_n}\right)}^2}& \left(38\right)\end{array}$

[0171] The received PDF model 154 is found by the convolution of the transmit PDF model with the noise PDF model

[0172] {tilde over (p)}_r(x)={tilde over (p)}_t(x)*p_n(x)  (39)

[0173] $\begin{array}{cc}{\tilde{p}}_r\left(x\right)={\int }_{-\infty }^{\infty }{\tilde{p}}_t\left(y\right)p_n\left(x-y\right)dy& \left(40\right)\\ {\tilde{p}}_r\left(x,{\sigma }_n\right)={\int }_{-\infty }^{\infty }\left\{({\mathrm{cW}}_{dc}\left(y\right)+\sum _{l=1}^L\frac{b_l}{\sqrt{2\pi }{\sigma }_l}{e}^{-{\left(\frac{x-m_l}{\sqrt{2}{\sigma }_l}\right)}^2})\left(\frac{1}{\sqrt{2\pi }{\sigma }_n}{e}^{-{\left(\frac{x-y}{\sqrt{2}{\sigma }_n}\right)}^2}\right)\right\}dy& \left(41\right)\end{array}$

[0174] where ${\sigma }_n^2=s/s_{\mathrm{nr}}=\mathrm{noise}\mathrm{variance}=\mathrm{noise}\mathrm{power}$

[0175] s=signal power

[0176] s_nr=signal-to-noise ratio (SNR)

[0177] {tilde over (p)}_r(x,σ_n)={tilde over (p)}_r(x,s_nr)=a function of both the received value and the noise power Carrying the SNR variable explicitly because it is the quantity of interest yields

[0178] {tilde over (p)}_r(x,s_nr)=α(x,s_nr)+β(x,s_nr)  (42)

[0179] where $\alpha \left(x,s_{\mathrm{nr}}\right)={\int }_{-\infty }^{\infty }\left\{({\mathrm{cW}}_{dc}\left(y\right))\left(\frac{1}{\sqrt{2\pi }{\sigma }_n}{e}^{-{\left(\frac{x-y}{\sqrt{2}{\sigma }_n}\right)}^2}\right)\right\}dy$ $\beta \left(x,s_{\mathrm{nr}}\right)={\int }_{-\infty }^{\infty }\left\{(\sum _{l=1}^L\frac{b_l}{\sqrt{2\pi }{\sigma }_l}{e}^{-{\left(\frac{y-m_j}{\sqrt{2}{\sigma }_l}\right)}^2})(\frac{1}{\sqrt{2\pi }{\sigma }_n}{e}^{-{\left(\frac{x-y}{\sqrt{2}{\sigma }_n}\right)}^2})\right\}dy$

[0180] The first term is integrated with a standard change of variables as $\begin{array}{cc}\alpha \left(x,s_{\mathrm{nr}}\right)=c\left[Q\left(\frac{-H_{\max }-x}{{\sigma }_n}\right)-Q\left(\frac{H_{\max }-x}{{\sigma }_n}\right)\right]{\sigma }_n^2=s/s_{\mathrm{nr}}& \left(43\right)\end{array}$

[0181] where [68, 70] $\begin{array}{cc}Q\left(\gamma \right)=\frac{1}{\sqrt{2\pi }}{\int }_{\gamma }^{\infty }{e}^{-\frac{1}{2}{\gamma }^2}d\gamma & \left(44\right)\end{array}$

[0182] An approximation to the Q-function has been found in that reduces the computation time by a factor of 2 over the MATLAB built-in error function (erf), which is related to the Q-function by a scale factor and change of variable. The approximation is $\begin{array}{cc}Q\left(x\right)\cong \left[\frac{1}{\left(1-\frac{1}{\pi }\right)x+\frac{1}{\pi }\sqrt{x^2+2\pi }}\right]\frac{1}{\sqrt{2\pi }}{e}^{-x^2/2}& \left(45\right)\end{array}$

[0183] and is accurate to within 1.5% across the entire range of x. Its use in the simulations helped control data run times.

[0184] The approximation has 1-1.5% accuracy for the argument in the range [0.3-1.2]. Closer to zero and greater than 1.2 yield more accurate approximations, asymptotically approaching zero error in both directions. For the windowed DC implementation here, the Q-function is used only in the first term of equation (42), where the positive argument shown in equation (43) results in the 1-1.5% error for the range $0.3<\frac{H_{\max }-x}{{\sigma }_n}<1.2H_{\max }-1.2{\sigma }_n<x<H_{\max }-0.3{\sigma }_n.$

[0185] Since Q(−x)=1−Q(x), the greatest inaccuracy is also in the range

[0186] 1.2σ_n−H_max>x>0.3σ_n−H_max.

[0187] The inaccuracy occurs in like manner for the other edge of the DC window. It is most pronounced just inside and outside the edges of DC window, and has greater contribution with larger channel noise power. The second term in the received PDF closed-form model equation (42) is exact, and reduces the contribution of the Q-function error to the received PDF convolution result.

[0188] The equation for the second term is solved using the result tabulated in a mathematical reference book, $\begin{array}{cc}{\int }_{-\infty }^{\infty }{e}^{\left(-p^2{x}^2\pm \mathrm{qx}\right)}dx=e^{\frac{q^2}{4p^2}}\frac{\sqrt{\pi }}{p}p>0& \left(46\right)\end{array}$

[0189] which results in $\begin{array}{cc}\beta \left(x,s_{\mathrm{nr}}\right)=\sum _{l=1}^L\frac{b_l}{\sqrt{2\pi }{\sigma }_{\mathrm{nl}}}{e}^{-{\left(\frac{x-m_l}{\sqrt{2}{\sigma }_{\mathrm{nl}}}\right)}^2}& \left(47\right)\end{array}$

[0190] where σ_nl²=σ_n²+σ_l²and ${\sigma }_n^2=s/s_{\mathrm{nr}}.$

[0191] The convolution of Gaussian functions results in a Gaussian function whose variance (power) is equal to the sum of the constituent variances (powers), which is a well-known result in probability theory.

[0192] FIGS. 19aand 19band 20aand 20billustrate the effectiveness of these equations in calculating PDF traces for 6 & 15 dB SNR and 30 & 60 dB SNR, respectively. In all cases, the top plot is an empirical data run of 1,000,000 iterates, the middle plot is the convolution result plotted from the equation

[0193] {tilde over (p)}_r(x,s_nr)=α(x,s_nr)+β(x,s_nr),  (48)

[0194] and the bottom plot is the superposition of the two.

[0195] Turning back to FIG. 18, the next step is to determine the maximum likelihood SNR 158. The basics of this method draw on the observation that the received probability {tilde over (p)}_r(x,s_nr) is a function of two variables, x and s_nr. It is not a joint probability function because only the function of x is a PDF with unity area (slices of constant s_nr), while the function of s_nrhas unbounded area (slices of constant x). FIGS. 19a, 19b, 20aand 20bshow the x dependency, which are PDFs.

[0196] Plots of the two-dimensional function are shown in FIGS. 21 and22, with slightly different angles of observation to highlight the variations in the functions with s_nrand enable good views of the peaks and valleys of the surface topology. These peaks and valleys became an important feature in the determination of the average SNR estimate. FIG. 23 illustrates examples of the PDF dependency on SNR for various values of the received iterate. These shapes are continuous and vary smoothly with SNR, which opened the possibility of an algorithmic method of solving for the maximum probability point to determine the corresponding SNR as the desired instantaneous value. The preferred method will be fast and respond only to maxima. In addition, the question of choosing local or global solutions arises in situations with more than one peak, as well as traces with two identical peaks that must be addressed.

[0197] Turning back to FIG. 18, either the Newton-Raphson method or the modified Newton-Raphson method is used to determine the maximum probability SNR 159. The Newton-Raphson iterative root approximation optimization technique was targeted for the algorithm because it converges quadratically to a solution, while all other possibilities such as bisection, secant, and regula-falsi converged linearly or slightly better than linearly. The Newton-Raphson technique finds the zeros of a function and converges to local solutions or diverges depending on its initial estimate and the correction step size. The individual iterate decisions are then smoothed through a running average calculation 160 and processing continues in FIG. 17.

[0198] The Newton-Raphson method operates on functions expressed as f(x_sol)=0, and requires the derivative of the function in its approximation correction process. Its iterated equation is $\begin{array}{cc}{x}_{j+1}=x_j-\frac{f\left(x_j\right)}{f^{\prime }\left(x_j\right)}& \left(49\right)\end{array}$

[0199] where

[0200] j=the iterate counter

[0201] x_j=the current approximation

[0202] f(x_j)=[f(x)]_x=x

[0203] $f^{\prime }\left(x_j\right)={\left[df\left(x\right)/dx\right]}_{x=x_j}$

[0204] x_j+1=the next approximation

[0205] At x=X_sol, f(x_sol)=0 and x_j+1=x_j=x_soland the iteration is complete.

[0206] Convergence of this method is not guaranteed and generally depends on the initial estimate and the step size that is calculated from the function and its derivative on each iterate. This method will find any zero of the function that is close enough to its initial value and has a well-behaved sequence of correction calculations. Although all references seen to date have similar descriptions of the technique that are not based on chaos theory, any solution point can be viewed from the standard nonlinear dynamical perspective as a fixed point of the iterative process. Correction steps that are large enough to move the next approximation out of the basin of attraction for a given solution will cause the iterative procedure to diverge from that point.

[0207] The function with zeros that the Newton-Raphson method solves for can be the derivative of another function, for example (cf. equation 42) $\begin{array}{cc}f\left(s_{\mathrm{nr}}\right)=\frac{d\left[{\tilde{p}}_r\left(x,s_{\mathrm{nr}}\right)\right]}{d{s}_{\mathrm{nr}}}& \left(50\right)\end{array}$

[0208] The solutions f(s_nr)=0 in this case correspond to the maxima, minima, and points of inflection for the probability function of SNR for a given received value, where its first derivative is zero.

[0209] Alternatively the modified Newton-Raphson technique may be used to determine the maximum probability SNR 159 by adapting the Newton-Raphson technique to this design via three main modifications to the control its functional characteristics. The first modification causes it to find only the maxima of the probability function of SNR based on the fact that the second derivative of a function identifies its concave down (maximum) or concave up (minimum) characteristic. The second change limits the maximum correction per step in order to avoid jumping over maxima without detecting them. And finally, the computed step size is multiplied by a factor that toggles between two values to prevent a period two fixed point of the iterative process spanning a peak and causing detection failure. Other standard programming techniques such as limiting the number of iterations allowed in each execution, hard programming upper and lower s_nrrange limits, protecting against a zero denominator, and early loop termination via testing each new approximation for proximity to the current value round out its implementation.

[0210] The Newton-Raphson equation expressed as functions of the probability function of SNR is $\begin{array}{cc}{s}_{{\mathrm{nr}}_{j+1}}=s_{{\mathrm{nr}}_j}-\frac{{\left[\frac{d\left[{\tilde{p}}_r\left(x,s_{\mathrm{nr}}\right)\right]}{d{s}_{\mathrm{nr}}}\right]}_{s_{{\mathrm{nr}}_j}}}{{\left[\frac{d{\hspace{0.17em}}^2\left[{\tilde{p}}_r\left(x,s_{\mathrm{nr}}\right)\right]}{d{s}_{\mathrm{nr}}^2}\right]}_{s_{{\mathrm{nr}}_j}}}& \left(51\right)\end{array}$

[0211] The derivative in the numerator is carried out below, and the fact that it equals zero at the solution allows a term to be dropped that simplifies the second derivative in the denominator. The correction term actually used is, therefore, not exactly as specified in equation (51), and so the notation will be changed to reflect this.

[0212] Let $\begin{array}{cc}{s}_{{\mathrm{nr}}_{j+1}}=s_{{\mathrm{nr}}_j}-\frac{f\left(s_{{\mathrm{nr}}_j}\right)}{f^{\prime }\left(s_{{\mathrm{nr}}_j}\right)}& \left(52\right)\end{array}$

[0213] where $f\left(s_{\mathrm{nr}}\right)\propto \frac{d\left[{\tilde{p}}_r\left(x,s_{\mathrm{nr}}\right)\right]}{d{s}_{\mathrm{nr}}}\mathrm{and}{f}^{\prime }\left(s_{\mathrm{nr}}\right)=\frac{df\left(s_{\mathrm{nr}}\right)}{d{s}_{\mathrm{nr}}}\propto \frac{d{\hspace{0.17em}}^2\left[{\tilde{p}}_r\left(x,s_{\mathrm{nr}}\right)\right]}{d{s}_{\mathrm{nr}}^2}$

[0214] Differentiating equation (42), $\begin{array}{cc}\frac{d\left[{\tilde{p}}_r\left(x,s_{\mathrm{nr}}\right)\right]}{d{s}_{\mathrm{nr}}}=\frac{d\left[\alpha \left(x,s_{\mathrm{nr}}\right)\right]}{d{s}_{\mathrm{nr}}}+\frac{d\left[\beta \left(x,s_{\mathrm{nr}}\right)\right]}{d{s}_{\mathrm{nr}}}& \left(53\right)\end{array}$

[0215] From equation (43) $\begin{array}{cc}\frac{d\alpha \left(x,s_{\mathrm{nr}}\right)}{d{s}_{\mathrm{nr}}}=c\left[\frac{d}{d{\sigma }_n}\left\{Q\left(\frac{-H_{\max }-x}{{\sigma }_n}\right)\right\}\frac{d{\sigma }_n}{d{s}_{\mathrm{nr}}}-\frac{d}{d{\sigma }_n}\left\{Q\left(\frac{-H_{\max }-x}{{\sigma }_n}\right)\right\}\frac{d{\sigma }_n}{d{s}_{\mathrm{nr}}}\right]& \left(54\right)\end{array}$

[0216] The differentiation of the Q-function was done using the Liebniz rule $\begin{array}{cc}\frac{d}{dy}\left[{\int }_{s\left(y\right)}^{r\left(y\right)}{p}_X\left(x\right)dx\right]=p_X\left(r\left(y\right)\right)\frac{dr\left(y\right)}{dy}-p_X\left(s\left(y\right)\right)\frac{ds\left(y\right)}{dy}.& \left(55\right)\end{array}$

[0217] The result was verified with the following equation $\begin{array}{cc}\frac{d}{dx}\mathrm{erf}\left(x\right)=\frac{2}{\sqrt{\pi }}{e}^{-x^2}& \left(56\right)\end{array}$

[0218] and the relationship $\begin{array}{cc}Q\left(x\right)=\frac{1}{2}\left[1-\mathrm{erf}\left(\frac{x}{\sqrt{2}}\right)\right]& \left(57\right)\end{array}$

[0219] as $\begin{array}{cc}\frac{d}{dx}Q\left(x\right)=-\frac{1}{\sqrt{2\pi }}{e}^{-\frac{x^2}{2}}& \left(58\right)\end{array}$

[0220] So $\begin{array}{cc}\frac{d\alpha \left(x,s_{\mathrm{nr}}\right)}{d{s}_{\mathrm{nr}}}=\left(\frac{1}{\sqrt{2\pi }{s}_{\mathrm{nr}}}\right)\left(\frac{c}{\sqrt{2}}\right)\left[\left(\frac{H_{\max }-x}{\sqrt{2}{\sigma }_n}\right)e^{-{\left(\frac{H_{\max }-x}{\sqrt{2}{\sigma }_n}\right)}^2}+\left(\frac{H_{\max }+x}{\sqrt{2}{\sigma }_n}\right)e^{-{\left(\frac{H_{\max }+x}{\sqrt{2}{\sigma }_n}\right)}^2}\right]& \left(59\right)\end{array}$

[0221] Carrying out the differentiation for the second term of equation (53), using equation (47) $\begin{array}{cc}\frac{d\beta \left(x,s_{\mathrm{nr}}\right)}{d{s}_{\mathrm{nr}}}=\sum _{l=1}^L\left\{\frac{d}{d{\sigma }_{\mathrm{nl}}}\left[\frac{b_l}{\sqrt{2\pi }{\sigma }_{\mathrm{nl}}}{e}^{-{\left(\frac{x-m_l}{\sqrt{2}{\sigma }_{\mathrm{nl}}}\right)}^2}\right]\left(\frac{d{\sigma }_{\mathrm{nl}}}{d{\sigma }_n}\right)\left(\frac{d{\sigma }_n}{d{s}_{\mathrm{nr}}}\right)\right\}& \left(60\right)\\ =\left(\frac{1}{\sqrt{2\pi }{s}_{\mathrm{nr}}}\right)\left(\frac{{\sigma }_n^2}{2}\right)\sum _{l=1}^L\left\{b_l{e}^{-{\left(\frac{x-m_l}{\sqrt{2}{\sigma }_{\mathrm{nl}}}\right)}^2}\left(\frac{1}{{\sigma }_{\mathrm{nl}}^3}\right)\left(1-{\left(\frac{x-m_l}{{\sigma }_{\mathrm{nl}}}\right)}^2\right)\right\}& \left(61\right)\end{array}$

[0222] Both terms of the differentiation for the Newton-Raphson correction term numerator have been found. Setting the numerator equal to zero allows the elimination of one s_nrterm, which simplifies the next differentiation for the denominator of the correction term. It is possible to form the total correction shown in equation (52) prior to simplification, but the leading term eliminated now because of numerator equality to zero contains the variable s_nr, which would enter into the differentiation process and result in more complicated equations. The notation change is shown below. $\begin{array}{cc}{\left[\frac{d\left[{\tilde{p}}_r\left(x,s_{\mathrm{nr}}\right)\right]}{d{s}_{\mathrm{nr}}}\right]}_{s_{{\mathrm{nr}}_{\mathrm{sol}}}}={\left[\frac{d\left[\alpha \left(x,s_{\mathrm{nr}}\right)\right]}{d{s}_{\mathrm{nr}}}\right]}_{s_{{\mathrm{nr}}_{\mathrm{sol}}}}+{\left[\frac{d\left[\beta \left(x,s_{\mathrm{nr}}\right)\right]}{d{s}_{\mathrm{nr}}}\right]}_{s_{{\mathrm{nr}}_{\mathrm{sol}}}}=0& \left(62\right)\end{array}$

[0223] Both of these terms have a common leading factor, which can be eliminated. $\begin{array}{cc}{\left[\frac{d\left[{\tilde{p}}_r\left(x,s_{\mathrm{nr}}\right)\right]}{d{s}_{\mathrm{nr}}}\right]}_{s_{{\mathrm{nr}}_{\mathrm{sol}}}}=\left(\frac{1}{\sqrt{2\pi }{s}_{{\mathrm{nr}}_{\mathrm{sol}}}}\right)\mathrm{f1}\left(s_{{\mathrm{nr}}_{\mathrm{sol}}}\right)+\left(\frac{1}{\sqrt{2\pi }{s}_{{\mathrm{nr}}_{\mathrm{sol}}}}\right)\mathrm{f2}\left(s_{{\mathrm{nr}}_{\mathrm{sol}}}\right)=0& \left(63\right)\end{array}$

[0224] Resulting in the new notation $\begin{array}{cc}f\left(s_{\mathrm{nr}}\right)=f_1\left(s_{\mathrm{nr}}\right)+f_2\left(s_{\mathrm{nr}}\right)=\sqrt{2\pi }{s}_{\mathrm{nr}}\left[\frac{d\left[{\tilde{p}}_r\left(x,s_{\mathrm{nr}}\right)\right]}{d{s}_{\mathrm{nr}}}\right]& \left(64\right)\end{array}$

[0225] where $\begin{array}{cc}{f}_1\left(s_{\mathrm{nr}}\right)=\left(\frac{c}{2}\right)\left[\left(\frac{H_{\max }-x}{{\sigma }_n}\right)e^{-\frac{1}{2}{\left(\frac{H_{\max }-x}{{\sigma }_n}\right)}^2}+\left(\frac{H_{\max }+x}{{\sigma }_n}\right)e^{-\frac{1}{2}{\left(\frac{H_{\max }+x}{{\sigma }_n}\right)}^2}\right]& \left(65\right)\\ f_2\left(s_{\mathrm{nr}}\right)=\left(\frac{{\sigma }_n^2}{2}\right)\sum _{l=1}^L\left\{b_l{e}^{-\frac{1}{2}{\left(\frac{x-m_l}{{\sigma }_{\mathrm{nl}}}\right)}^2}\left(\frac{1}{{\sigma }_{\mathrm{nl}}^3}\right)\left(1-{\left(\frac{x-m_l}{{\sigma }_{\mathrm{nl}}}\right)}^2\right)\right\}& \left(66\right)\end{array}$

[0226] The denominator in equation (52) becomes

[0227] f′(s_nr)=f₁′(s_nr)+f₂′(s_nr)  (67)

[0228] with $\begin{array}{cc}{f}_1^{\prime }\left(s_{\mathrm{nr}}\right)=\left(\frac{1}{2s_{\mathrm{nr}}}\right)\left(\frac{c}{2}\right)\left[\left(\frac{H_{\max }-x}{{\sigma }_n}\right)e^{-\frac{1}{2}{\left(\frac{H_{\max }-x}{{\sigma }_n}\right)}^2}\left\{1-{\left(\frac{H_{\max }-x}{{\sigma }_n}\right)}^2\right\}+\left(\frac{H_{\max }+x}{{\sigma }_n}\right)e^{-\frac{1}{2}{\left(\frac{H_{\max }+x}{{\sigma }_n}\right)}^2}\left\{1-{\left(\frac{H_{\max }+x}{{\sigma }_n}\right)}^2\right\}\right]& \left(68\right)\\ f_2^{\prime }\left(s_{\mathrm{nr}}\right)=\left(\frac{{\sigma }_n^2}{2s_{\mathrm{nr}}}\right)\left(\frac{{\sigma }_n^2}{2}\right)\sum _{l=1}^L{b}_l{e}^{-\frac{1}{2}{\left(\frac{x-m_l}{{\sigma }_{\mathrm{nl}}}\right)}^2}\left(\frac{1}{{\sigma }_{\mathrm{nl}}^3}\right)\left\{\left(\frac{1}{{\sigma }_{\mathrm{nl}}^2}\right)\left[3{\left(1-{\left(\frac{x-m_l}{{\sigma }_{\mathrm{nl}}}\right)}^2\right)}^2-2{\left({\left(\frac{x-m_l}{{\sigma }_{\mathrm{nl}}}\right)}^2\right)}^2\right]-\frac{2}{{\sigma }_n^2}\left(1-{\left(\frac{x-m_l}{{\sigma }_{\mathrm{nl}}}\right)}^2\right)\right\}& \left(69\right)\end{array}$

[0229] The equations are now fully defined, and the algorithmic modifications can be made. Equation (52) is repeated below as the starting point, and depicts the Newton-Raphson method with attractive fixed points at all maxima, minima, and points of inflection for the probability function {tilde over (p)}_r(x,s_nr). $\begin{array}{cc}{s}_{{\mathrm{nr}}_{j+1}}=s_{{\mathrm{nr}}_j}-\frac{f\left(s_{{\mathrm{nr}}_j}\right)}{f^{\prime }\left(s_{{\mathrm{nr}}_j}\right)}& \left(70\right)\end{array}$

[0230] The minima have a concave up structure with a positive second derivative of {tilde over (p)}_r(x,s_nr), while the maxima have a concave down structure and a negative second derivative. The same holds true for the corresponding portions of a point of inflection in that the side leading up into the point is concave down and the side leading up away from the point is concave up. A simple test prior to the execution of equation (70) turns the minima into repelling fixed points, while maintaining the maxima as attractive fixed points.

[0231] If f′(s_nr)>0% Concave up−minima

[0232] f′(s_nr)=−f′(s_nr) % Negates correction direction  (71)

[0233] % therefore pushes away from minima

[0234] end $s_{{\mathrm{nr}}_{j+1}}=s_{{\mathrm{nr}}_j}-\frac{f\left(s_{{\mathrm{nr}}_j}\right)}{f^{\prime }\left(s_{{\mathrm{nr}}_j}\right)}$

[0235] The maximum probability algorithm performs better when working with SNR values expressed as dB rather than ratios. Although the typical calculation method is to work with power ratios and convert results to dB, the traces shown in FIG. 23 have an intuitively good appearance similar to polynomial curves on a dB SNR scale, so the various considerations are expressed in dB for this situation.

[0236] A set of maximum step criteria were established to avoid jumping over maxima without detecting them. A sliding maximum step was used because the areas of these traces that showed little or no activity could tolerate larger steps than those areas with a lot of peaks and valleys, and maximized the efficiency of the algorithm at locating the peaks. The steps and their regions of definition are shown in Table 3.

[0237] The maximum step size is multiplied at each pass by a variable that toggles between the values of 1 and ⁵½=0.8706, which is an empirically determined value. The roots of ½ are used as candidates for a suitable step size modification because the various roots yield conveniently spaced possibilities. A number of other selection methods would have worked as well. The dual step size technique avoided the situation of a period two fixed point causing the algorithm to perpetually span the desired maximum point without ever finding it within the SNR range of a given maximum step size. Care was required at the crossover points of the sliding maximum steps that were used, where the wrong choice of toggle multiplier could map the larger range exactly into the smaller and cause algorithm failure due to a period two fixed point.

[0238] The detection rails were set at 0 dB and 100 dB SNR, which were the expected minimum SNR of performance and essentially no noise, respectively. Next approximation values that exceeded these limits were adjusted, so instantaneous SNR values are always within this range. These were chosen prior to completion of the receiver design and the subsequent evaluation of its performance, and especially the 0 dB minimum SNR may be a good subject for further investigation.

[0239] At the completion of the Newton-Raphson computation, the instantaneous maximum probability SNR estimate for iterate k has been found.

[0240] s∃_nr[k]=s_nr=[s_nr]_{s−s<ε} where ε=desired accuracy  (72)

[0241] Turning back to FIG. 17, the next step in the receiver estimation engine processing is to calculate the running weighted average 115.

[0242] Combining the instantaneous SNR estimates on every received iterate into an average value that accurately tracks is accomplished by a running weighted average calculation whose weights are derived from probability computations so the estimated SNR tracks the true SNR closely. The current received iterate and the corresponding instantaneous SNR estimate are processed to find their probability of occurrence. The bin width w is ignored because it is the same for all bins in the PDF model and contributes nothing to the decision process. Only the bin values of the PDF functions will be used in the weighted averaging process and will be loosely referred to as the probability. For each received iterate, the SNR is estimated, and the probability of occurrence is found for the (received value, SNR estimate) pair. These values correspond to a point in the two-dimensional probability function shown in FIGS. 21 and 22. Both the received value and the SNR estimate will change for every iterate, so the instantaneous probability moves with respect to both axes on every calculation.

[0243] In order for the estimated SNR to track the true SNR, the averaging weights are defined as the inverse instantaneous probability of the (received value, SNR estimate) pair $\begin{array}{cc}{w}_{\mathrm{av}}\left[k\right]=1/P_i\left[k\right].& \left(74\right)\end{array}$

[0244] Two considerations for calculating the weighting factors are avoiding division by zero and limiting the maximum weight. Observed probabilities of zero within the resolution of MATLAB did happen, and the probability result is now restricted to a minimum value of 10⁻¹⁰to avoid division by zero. The maximum weight calculation limits the amount of emphasis put on the very low probability values so that their inverses do not swamp the average, and was derived from the observed average probability and the averaging window length.

[0245] {overscore (P)}_i=0.4444 $\begin{array}{cc}{w}_{\max }=\frac{1}{{\stackrel{\_}{P}}_i}\frac{W_{\mathrm{ra}}}{f_W}& \left(75\right)\end{array}$

[0246] where W_ra=the running average window length

[0247] and f_W=3, an empirical factor.

[0248] Data was taken for various combinations of variables, and the empirical factor value of three was found to yield good results. Averaging window lengths of 8 and 16 are responsive to changes in channel conditions but have a grassy appearance with some areas of significant discrepancy between the estimated SNR and the true SNR. Windows of 32 and 64 iterates showed smoother traces with better estimate tracking. This simulation was done with a 64-iterate window, and examples of its performance are shown in FIG. 24. In FIG. 24, it can be seen that the running weighted average SNR locks onto a close estimate of the true SNR within approximately 200 iterates for a 64-iterate window. The traces have some variation because they are estimates and not power measurements, however they maintain a close proximity to the true value once lock is attained. In simulations of the receiver, the receiver decisions are ignored until the 200-iterate transient period has passed, after which time data bit recovery commences. In the event that the channel SNR varies, the SNR estimator will track the changes in the same manner as the receiver initialization process did with an unknown SNR. Sharp variations will cause a 200-iterate transient period to occur, as happens at receiver initialization, and gradual changes will be tracked on a continual basis.

[0249] The issues of local versus global instantaneous SNR maximum likelihood detection and algorithm initialization are related. Global detection is desired in the absence of knowledge at estimate algorithm inception, but the average value drives the instantaneous peak detection once the averaging window has filled up with globally detected maximum likelihood instantaneous SNR values. This feedback of the current average SNR estimate serves as the initial approximation in the Newton-Raphson root approximation optimization procedure. The average SNR feedback enables the determination of local maxima in the vicinity of the current average SNR estimate after the its value has been established via global values, which tends to drive the estimated SNR toward an accurate assessment of the true SNR. Observations of the probability vs. SNR functions led to the empirical initial SNR value of 23 dB for the instantaneous SNR detector Newton-Raphson algorithm during the averaging window initialization. This value is located about midway between the lobes of the observed double-peak functions and had good mobility toward either higher or lower values for the various observed traces.

[0250] The techniques used in this SNR estimator do not require differentiable equations for the probability function to be successful implementations. Standard derivative definition equations could be used in the Newton-Raphson iteration technique where two close values are used with the corresponding function responses, and the second derivative can be found in the same manner. $\begin{array}{cc}\frac{df\left(x\right)}{dx}\cong \frac{f\left(x_2\right)-f\left(x_1\right)}{x_2-x_1}\mathrm{and}\frac{d{\hspace{0.17em}}^{2}f\left(x\right)}{d{x}^2}\cong \frac{\frac{df\left(x_2\right)}{dx}-\frac{df\left(x_1\right)}{dx}}{x_2-x_1}& \left(76\right)\end{array}$

[0251] It is necessary in the averaging process to find the probability of occurrence for the (received iterate, maximum probability SNR) pair, which can be done with empirical data rather than via closed-form equations. The difficulty with this approach is that the probability response as a function of SNR is not known and cannot be verified to be smooth and well behaved for proper Newton-Raphson iteration performance.

[0252] Turning back to FIG. 17, the transmit value Maximum a Posteriori (MAP) estimator 117 generates the second estimate for a given received iterate, which is the optimal probabilistic detector because it minimizes the probability of error. Differences between conventional communication methods and this chaotic system in the transmission characteristics, the modulation technique, and the presence of chaotic dynamics require additional considerations in the decision process. The MAP decision resulting in the second estimate is a statistical method that yields the lowest bit error rate.

[0253] The chaotic communication system herein relates the binary possibilities of logical 1 and logical 0 with a constantly changing sequence drawn from a continuum of possible values, the range of which is identical for both logical states. There is not a direct relationship between a transmit value and its message content. Any given transmit value could represent either logical state depending on the attractor from which it was generated at that point in time, and the next time the value (or a value within a small neighborhood) is encountered it may have come from the opposite attractor.

[0254] The MAP calculation is performed on every received iterate, so the iterate index k is carried explicitly in the following equations. In order to relate the continuum of values x_tto the PDF as empirically determined or modeled, reference is made to the formulations derived previously. There the empirical PDF traces were generated by collecting the countably infinite set of transmit values into quantized data bins of acceptably narrow widths. The PDF traces were then modeled by a set of weighted Gaussian functions, resulting in equation (25). This quantized data bin approach is extended below to the maximum a posteriori (MAP) statistical calculation, and a decision criterion is derived. This criterion is shown to be computationally intensive, and a simplified calculation that minimizes the computational load for a MAP calculation is developed to enable faster execution of the receiver estimation engine decision process.

[0255] The continuous range of countably infinite values is considered instead to have a very large but finite set l of possible values, as is encountered by a chaotic map processed on a computer because of the finite resolution of computer-generated numbers. The general relationship between the probability and a specific value in the range of a PDF is $\begin{array}{cc}P\left[A<X\le B\right]={\int }_A^B{p}_X\left(x_0\right)d{x}_0& \left(77\right)\end{array}$

[0256] The probability that the l^thtransmitted value possibility x_t(l) falls into bin b in the quantized PDF of section 4.2.2 is equal to $\begin{array}{cc}{P}_t\left[\left(x_{\mathrm{tb}}-w/2\right)<x_t\left(l\right)\le \left(x_{\mathrm{tb}}+w/2\right)\right]={\int }_{x_{\mathrm{tb}}-w/2}^{x_{\mathrm{tb}}+w/2}{}_{b}p_t\left(x_t\left[l\right]\right)d{x}_t=w\left[{}_{b}p_t\left(x_t\left[l\right]\right)\right]& \left(78\right)\end{array}$

[0257] where x_*b=* iterate series PDF central value of bin b (i.e.—t→Tx series,

[0258] r→Rx series, etc.)

[0259] _bp_*(●)=the probability of the bin into which the value (●) falls, for the * iterate series

[0260] l=a dummy variable of delineating the range of possible transmit values

[0261] and

[0262] b=1,2, . . . , B bin index.

[0263] Note that the index “l” used here is an indicator of the countably infinite possible transmit values, the total set of which numbered about 2.6×10¹⁴on the hardware platform used with the software package MATLAB. This set of transmit values became quantized into the 4096 bins selected as adequate to represent the discretized PDF with sufficient resolution. It is distinct from the index “l” used previously to enumerate the weighted Gaussian functions comprising the PDF models. In shorthand notation,

[0264] P_t(x_tb[l])=w[_bp_t(x_t[l])]  (79)

[0265] The determination of a MAP transmit value seeks the maximum value of the equation $\begin{array}{cc}\underset{b}{\max }\left\{P\left(x_{\mathrm{tb}}\left[l\right]x_r\left[k\right]\right)\right\}=\underset{b}{\max }\left\{{}_{b}p_{r|t}\left(x_r\left[k\right]|x_{\mathrm{tb}}\left[l\right]\right)P_t\left(x_{\mathrm{tb}}\left[l\right]\right)\right\}& \left(80\right)\\ =\underset{b}{\max }\left\{{}_{b}p_{r|t}\left(x_r\left[k\right]|x_{\mathrm{tb}}\left[l\right]\right)w\left[{}_{b}p_t\left(x_t\left[l\right]\right)\right]\right\}& \left(81\right)\end{array}$

[0266] where k is the time index and x_r[k] is the k^threceived value. The bin width term w is the same constant value for all probability bins, and so can be dropped from the maximization process. The resulting equation to maximize is $\begin{array}{cc}\underset{h}{\max }\left\{P\left(x_{\mathrm{th}}\left[l\right]|x_r\left[k\right]\right)\right\}=\underset{b}{\max }{\{}_b{p_{r|t}\left(x_r\left[k\right]|x_{\mathrm{th}}\left[l\right]\right)}_h{p}_t\left(x_t\left[l\right]\right)\}& \left(82\right)\end{array}$

[0267] The maximization occurs across all bins in the quantized PDF range of values. The index b can be dropped, with the understanding that the entire PDF is included in the maximization range, and that the resulting values will be quantized. A further notation simplification can be made by recognizing that the transmit PDF is fixed, regardless of which time iterate index k is current. In addition, the countably infinite number of possible transmit values can be considered continuous, allowing elimination of the index “l”. The results are

[0268] p_t(x_t[l])→p_t(x_t) is the transmit PDF  (83)

[0269] x_tb[l]→x_tthe l^thtransmit value possibility that  (84)

[0270] results in x_r[k] is the entire transmit range $\begin{array}{cc}\underset{x}{\max }\left\{P_{t|r}\left(x_t|x_r\left[k\right]\right)\right\}=\underset{x}{\max }\left\{p_{r|t}\left(x_r\left[k\right]|x_t\right)p_t\left(x_t\right)\right\}.& \left(85\right)\end{array}$

[0271] Equation (85) is the MAP technique for the estimation of continuous variables. The term p_r|t(x_r[k]|x_t) is a window that multiplies the transmit PDF p_t(x_t), and the maximization of the result is the desired MAP transmit value for the k^thsystem iterate. This window serves to define the region of the PDF over which the maximization will occur, and also to modify the shape of the transmit PDF according to the noise contained in the received sequence x_r.

[0272] FIG. 25 shows a flowchart of the general MAP transmit value determination. The transmit sequence PDF is determined and stored by empirical observation or using model calculations at receiver initialization 250. The channel noise characteristics for every received iterate are determined and the channel noise PDF is obtained by empirical observation or model calculations 251. A basic window is constructed with the channel noise PDF reversed along the random variable axis 252. The iterate widow is constructed by centering the basic window on the currently observed received value 253. The transmit PDF is multiplied by the iterate window 254. The maximum PDF value in the windowed transmit PDF is then found 255. The MAP estimate is the random variable value corresponding to the maximum PDF value in the windowed transmit PDF 256.

[0273] For the current situation with Gaussian noise and a Gaussian approximation to the transmit PDF, it is conceptually possible to perform this multiplication and get a closed-form equation for the result. This expression can then be differentiated once or twice, depending on the iterative root approximation modification algorithm used to find the maxima, as is done for the SNR estimator. However, the multiplication result will have many sharply defined peaks for the noise levels of interest for practical communications, and this approach was discarded as impractical. The multiplication result is more easily searched for its largest value with a software routine that locates the maximum value of a vector of numbers, which yields the most likely transmitted value. A ‘max’ function is standard in many mathematical platforms such as MATLAB, and makes the search for the maximum value in a vector of numbers routine.

[0274] The transmit PDF window construction is outlined below. The general result is that the noise PDF, reversed in x, is centered on the received value. This window function then multiplies the transmit PDF, and the maximum value of the result is found. A simplified decision calculation has been developed for this design, to ease the computational burden of the technique by avoiding the repeated computation of Gaussian exponentials for the determination of the MAP transmitted value.

[0275] Turning back to FIG. 17, the next step is to calculate the transmit PDF window 118. The fact that a continuum of values can be transmitted for either logical binary state complicates the determination of an a posteriori maximum likelihood transmit value, compared to the methods used in the prior art. Furthermore, the result of such a calculation is not necessarily the best estimate of the transmit value because the effects of chaotic dynamics have not been taken into account. In this implementation the MAP computation uses only the current value and makes no use of past or future received values. However, it is a good contributor to the receiver decision of the transmit value, along with some other estimates that help track the deterministic random motion of the chaotic system.

[0276] A graphical proof of the transmit PDF window construction is presented in FIG. 28. It draws on the fact that for additive, statistically independent noise, the received PDF is the convolution of the transmit PDF with the channel noise PDF

[0277] p_r(x)=p_t(x)*p_n(x).  (86)