Constrained pre­equalization accounting for multi­path fading emulated using large RC networks: applications to wireless and photonics communications equalization accounting multi­path fading emulated using large

Multi-path propagation is modelled assuming a multi-layer RC network with randomly allocated resistors and capacitors to represent the transmission medium. Due to frequency-selective attenuation, the waveforms associated with each propagation path incur path-dependent distortion. A pre-equalization procedure that takes into account the capabilities of the transmission source as well as the transmission properties of the medium is developed. The problem is cast within a Mixed Integer Linear Programming optimization framework that uses the developed nominal RC network model, with the excitation waveform customized to optimize signal ﬁdelity from the transmitter to the receiver. The objective is to match a Gaussian pulse input accounting for frequency regions where there would be pronounced fading. Simulations are carried out with different network realizations in order to evaluate the sensitivity of the solution with respect to changes in the transmission medium mimicking the multi-path propagation. The proposed approach is of relevance where equalization techniques are difﬁcult to implement. Applications are discussed within the context of emergent communication modalities across the EM spectrum such as light percolation as well as emergent indoor communications assuming various modulation protocols or UWB schemes as well as within the context of space division multiplexing.


Introduction
Multi-path propagation is a phenomenon that is often encountered when matter is probed by electrical acoustic or electromagnetic signals with wavelengths smaller than the physical dimensions of the features that need to be characterized (Lowery et al. 2012;Kaatze 2008). Interpretation of measurements where percolation of the excitatory wave through the material has taken place is never straightforward and requires adoption of solutions developed by the inverse problems community. It is also a typical problem encountered in communications where performance of wireless and mobile communication links is compromised by multi-path interference and frequency selective fading (Taub and Schilling 1986;Sari et al. 1995). The phenomenon is most prominent in indoor communication systems based on visible or IR optoelectronics based transceivers (Hashemi 1993) and limit proliferation of communications solutions at 60 GHz or THz links (Smulders 2002;Piesiewicz et al. 2005;Galvão et al. 2007a, b). In addition, it is encountered in optical fibre or waveguide based communications under multimode propagation and when photonic crystal fibres are used. It is generally accepted that irrespective of the propagation medium, most commonly adopted communication protocols such as orthogonal frequencydivision multiplexing (OFDM) employing multilevel modulation schemes with non-constant amplitude (e.g. 16 QAM) generally require estimation and tracking of the fading channel parameters to perform coherent demodulation (Li et al. 1998;Edfors et al. 1998;Yang et al. 2001). For particular modulation modalities, channel estimation is most commonly achieved by exploiting the correlation of the channel frequency response at different frequencies and times. Alternatively, channel estimation and multi-path fading may be based on a parametric modeling approach where the channel is modeled as a finite impulse response (FIR) filter. A channel estimator in its broader sense is often used to estimate the time delays, gains, and phases of the paths.
In the following sections, we first discuss the relevance of three-dimensional models of RC input-output network responses to communications as this provides a novel way of modelling multi-path propagation. In the proposed approach, multi-path propagation is seen as analogous to a percolation process. We further explain that frequency selective attenuation can be parsimoniously described using system identification models which can be associated with the proposed network structure. This approach is based on the physical insight that has been recently developed by studying the underlying percolation processes discussed in the Physics literature (Stroud 1979;Soukoulis et al. 1987;Ho et al. 1990;Rivas et al. 1999;Hu et al. 2004).
Because of the universal relations underlining coupling of optical power between microresonators and dielectric waveguides (Yariv 2000), and the fact that interaction of multiconductor transmission lines with electromagnetic radiation is usually accounted for through distributed current and voltage sources along the lines (with the current sources being related to the electric field and the voltage sources being proportional to the magnetic field), the use of percolation models for electrical networks (Desoer and Kuh 1969) should also be seen as complementing modal transmission line theory (Snyder and Love 1983;Hoefer 1985;Cangellaris 1988;Faria 1993;Castellanos et al. 1997;Rowe 1999;Akbari et al. 2000;Lin et al. 2002;Pregla 2008;Nitsch et al. 2009;Stutzman 2012;Wait 2013;Rambousky et al. 2013;Spadacini et al. 2014). As discussed in a recent review (Kafesaki and Soukoulis 2006), the propagation of EM waves after taking into consideration the transmission properties of the medium can be theoretically studied using the plane wave (PW) method, the transfer matrix method (TMM) or the finite difference time domain method (FDTD). Such simulations are normally complemented by using mode solvers such as CUDOS MOF Utilities (http://sydney.edu.au), RSoft (http://optics. synopsys.com) or COMSOL (http://www.comsol.com). The above background provides a justification for the adoption of the proposed networks to model multi-path propagation.
An advantage of the proposed approach is that it enables the development of equalization techniques that would be of relevance across the communications industry (radio, mm-wave, terahertz or photonic based), assuming coherent, incoherent or partially coherent emitter modalities and homodyne or heterodyne receiver modalities. Channel dynamics for multiple paths are written in a state space form which can be converted into an input-output realization where Laplace transforms can be used to describe frequency selective attenuation in a parsimonious manner. We then develop a pre-equalization procedure where the objective is to match the output waveform with the input function. In our case a Gaussian function is used to achieve this. The difficulty in the task lies on the fact that the proposed approach should account for constraints on instantaneous values of voltage, current and power that can be either delivered by the source or absorbed by the transmission medium. The problem is cast within a Mixed Integer Linear Programming (MILP) optimization framework that uses the developed nominal RC network model, with the excitation waveform customized to optimize signal fidelity from the transmitter to the receiver. Simulation results are then provided, showing the robustness of the proposed formulation. A general discussion of other applications in photonics and measurement science that would benefit from the proposed formulation is also provided.

Channel modelling using RC networks
In this paper, we model multi-path propagation through a transmission medium by using large three-dimensional (3D) RC networks. The rationale behind this type of model is based on the fact that for a simple RC network, with voltage input uðtÞ and current iðtÞ flowing through a resistor R and a capacitor C, the dynamic relation between the input voltage and the current can be described by a transfer function of the form: where =[.] denotes the Laplace transform.
If the signals uðtÞ and iðtÞ are measured at sampling times kT, where T is a sampling period and k is the sample index, the relation between sequences uðkTÞ and iðkTÞ can be expressed by a transfer function of the form: where Z[.] denotes the Z-transform. The values of a and b depend on the sampling period T, as well as the discretization method employed. If the input voltage is determined by a digital processor, then its value can be expected to remain constant over the sampling period associated to the operation of the digital-to-analog converter. In such a case, the zero-order-hold (ZOH) method can be employed (Dickinson 1991), which results in a ¼ exp½ÀTðRCÞ À1 . The value of b can be obtained by noting that the ZOH method preserves the DC gain of the transfer function, i.e. GðsÞj s¼0 ¼ G d ðzÞj z¼1 . Since GðsÞj s¼0 ¼ 0 in this case, it follows that b ¼ 1 À a ¼ 1 À exp½ÀTðRCÞ À1 . It can easily be seen that 0\a\1 and 0\b\1. Finally, the transfer function G d ðzÞ can be expanded in a power series of z as: It is worth noting that the term z Ài in the Z-transform domain is associated to a delay of i sampling periods in the time domain. Therefore, the dynamics of the network can be interpreted as a combination of a direct (line of sight) transmission path and multiple reflections with gains that decrease as the associated delay is increased (Reis and Galvão 2004). Such an interpretation is consistent with multipath propagation, as the attenuation tends to be stronger over longer paths.
A large 3D RC network such as the one shown in Fig. 1 is employed to represent a complex medium, comprising a multitude of propagation paths with different delays and attenuations. A typical approach to handle the distortions of the transmitted signal caused by multi-path propagation is to include a channel equalizer at the receiver side. In the present case, if uðtÞ and iðtÞ are regarded as the transmitted and received signals, respectively, a discrete-time equalizer could be implemented in the form of a finite impulse response (FIR) filter as: with weights w 0 ; w 1 ; . . .w L that should be adjusted according to the dynamics of the transmission medium. Such an adjustment can be carried out in an adaptive manner by using the well-known Least-Mean-Squares (LMS) algorithm e.g. (Widrow and Stearns 1985). The LMS adaptation rule can be expressed as: where l [ 0 is an adaptation parameter, wðkÞ and iðkTÞ denote the vectors of weights and inputs of the FIR filter and at kth sampling time: and eðkTÞ denotes the difference between the filter output and the desired response r(kT): It is worth noting that this adaptation scheme requires the knowledge of the actual transmitted signal. In practice, the adaptation can be carried out by transmitting a known ''training'' signal. Alternatively, blind equalization methods may be employed (Johnson et al. 1998).
The present paper explores a different approach, where the shape of the transmitted signal is optimized by taking into account the dynamics of the medium and the capabilities of the transmitter source, in order to achieve a desired waveform at the receiver. The proposed approach can be regarded as a pre-equalization procedure, which dispenses with the need for the transmission of training signals. This has applications in RF wireless (Koenig et al. 2013) or light communications (Chow et al. 2012;Ghassemlooy et al. 2013) or in coherent optical receivers for optical fibres based communications (Savory 2010;Mori et al. 2013;Chen and Chng 2005;Liu et al. 2014). Furthermore, since such networks can be emulated in real-time using reconfigurable FPGAs, there are likely to be additional applications in FPGA based digital coherent optical OFDM signal processing which offer new opportunities for spectral efficiency, receiver sensitivity, and polarization-dispersion resilience (Yang et al. 2009).
Furthermore, it is worth noting that the proposed multi-path models underlying the proposed method should be particularly useful in developing equalization for spatial mode based multiplexing schemes which are likely to provide additional capacity in future Fig. 1 Schematic diagram of a 3D RC network with random allocation of R, C elements. The grey plates indicate a pair of electrodes employed to connect the network to a voltage source. The output resistance of the source is denoted by R S generation communication modalities (Richardson et al. 2013;van Uden et al. 2014;Mizuno et al. 2014;Huang et al. 2015;Shieh et al. 2012). Because space division multiplexing is also likely to be important in astronomical applications (Leon-Saval et al. 2010) as well as on-chip high throughput optical communications (Luo et al. 2014), the proposed modelling procedure accounting for multi-path propagation with the excitation waveform customized to optimize signal fidelity from the transmitter to the receiver should have a wide domain of technological applications.

Optimization of the input waveform
The problem addressed herein consists of determining an appropriate waveform for the input voltage uðtÞ in order to obtain a current waveform iðtÞ that is closest to a desired profile, respecting constraints on instantaneous values of voltage, current and power that can be either delivered by the source or absorbed by the transmission medium. More formally, we wish to minimize a 1-norm cost function of the form subject to: where rðkTÞ; k ¼ 0; 1; . . .; N, is the desired current profile, specified over N sampling times, and u max , i max , P max are bounds on the voltage, current and power of the source. It is assumed that the RC network is initially at rest (uncharged capacitors) and that the relation between the sequences uðkTÞ and iðkTÞ; k ¼ 0; 1; . . .; N is expressed by a model derived from the topology and R, C component values of the network. More specifically, we have adopted the modelling framework presented in  to derive a linear, continuous-time state-space model, which was subsequently simplified by using an orderreduction procedure based on balanced realizations. A review of model order reduction techniques that would be suitable for the task is found in (Heydari and Pedram 2006). For use in the optimization problem, the reduced-order model was then discretized in the form: where xðkTÞ 2 < n denotes the n-dimensional state vector at the kth sampling time, A, B, C are matrices of dimensions (n9n), (n91), (19n), respectively, and D is a scalar. For more details concerning the discretization of state-space models, the reader is referred to (Franklin et al. 1997).
The minimization of the 1-norm cost function in (8) subject to the dynamic relations expressed by the model Eq. (10a, b) and the bounds on the voltage (9a) and current (9b) values can be cast into the form of linear programming (Camacho and Bordons 1999). However, the presence of the power constraints (9c) result in a non-convex set of feasible ðuðkTÞ; iðkTÞÞ values. In the present work, this issue is handled by using linear approximations that result in a mixed-integer linear programming (MILP) problem, as detailed below. In what follows, the indication of the sampling period T will be omitted in the equations for brevity.

Conversion of the optimization problem to a MILP problem
As stated in the previous section, the inclusion of power constraints entails an optimization problem over a non-convex region of possible solutions. This issue poses a difficulty to obtain the optimal solution, as most optimization algorithms are designed to search over a convex set of feasible solutions. Therefore, it is necessary to restate the optimization problem to an equivalent form, which in turn can be solved with standard convex optimization techniques. Among the alternatives to recast the optimization problem is the conversion to a Mixed Integer Programming (MIP) problem involving both continuous and binary-valued variables.
A procedure that can be employed to convert the ''or''-type constraints that result in a non-convex optimization problem to the ''and''-type which entail convex optimization with the presence of binary variables, is the ''big-M'' method (Richards et al. 2002). A conversion of this type was presented in the formulation by Lima et al. 2010, which treated the case of power constraints in an electrical motor system by approximating the power constraints by straight line segments. In that work, the cost function was quadratic in terms of the optimization variables thus resulting in a Mixed Integer Quadratic Programming (MIQP) problem. Contrary to that work, in the optimization required in the current study, the cost function is linear in terms of the optimization variables, and should be casted as a Mixed Integer Linear Programming (MILP) problem. MILP problems are less intensive in terms of computational effort than MIQP ones, therefore, if one can obtain a linear cost whose optimization reflects the desired behavior of the system in a satisfactory manner, the computational effort to obtain the optimal solution may be reduced.
MIP problems present the drawback of exponential growth of the computational complexity with the increase of the number of binary variables. As a consequence, in obtaining a problem that can be solved in a reasonable amount of time as would be required in real-time pre-equalization applications, it is important to keep the number of binary decision variables to a minimum. In (Lima et al. 2010) the ''or''-type constraints were written in terms of the bounds of the feasible region, therefore, one constraint corresponded to each half-plane of the feasible region. In contrast, in the present work the constraints are expressed in terms of convex polytopes whose union is the feasible region. In general, the number of such polytopes, N p , is less than the number of half-planes, allowing for a reduction in the number of ''or''-type constraints. Figure 2 illustrates the benefits of using the polytopes for the constraints as opposed to the half-planes. As can be seen, there are 22 half-planes that need to be used with ''or''-type constraints to represent the feasible region, in contrast to 21 convex polytopes. In this figure, only the power constraints are depicted, as they require an approximation by polytpes. The bounds on the voltage and current are not shown, given that their representation is a simple rectangle parallel to the axes.
In the traditional employment of the ''big-M'' method, each ''or''-type constraint corresponds to one binary variable in the optimization problem. In a recent paper (Prodan et al. 2012), the authors presented an extension to that method to assign tuples of binary variables to each constraint, allowing for a number of binary variables N b ¼ log 2 N p . By implementing this tuple association to feasible regions in this work, the number of binary variables necessary to represent the non-convex constraints is considerably reduced. The application of such method, in turn, entails a number of additional constraints on the binary variables to render the possibly unallocated tuples infeasible. This number of additional constraints is kept to a minimum by following the procedure developed in (Afonso and Galvão 2014). Therefore, the constraints over the current and voltage can be written as: In (11), the N p convex polytopes C j ; 1 j N p , are described by the pairs of matrices and vectors S j and h j , respectively. F j is the number of facets of the polytope j and 1 Á Á Á 1 ½ T Fj is a vector of F j 9 1. Each function f j is zero only if the tuple k corresponds to the value that was assigned to render the constraints active, otherwise f j is equal to a positive integer. Since the constant scalar M is large enough to render all inequalities inactive, only one polytope constraint is enforced at each time step, as the others are relaxed: where L is a positive integer. To render the unassigned tuples infeasible, one constraint is appended to the optimization problem for every time instant k: where b l is the l-th binary variable of the tuple, in accordance to Afonso and Galvão (2014). The cost function that was adopted in this work is the 1-norm of the vector formed by stacking the difference between the desired current value and the predicted one at every sample time, according to Eq. (8). To achieve this, auxiliary continuous variables are employed as means to penalize the absolute values of the deviations: by imposing the following additional constraints: nðkÞ ! rðkÞ À iðkÞ which ensure that nðkÞ ¼ jrðkÞ À iðkÞj upon minimization of (14). The resulting optimization problem is: subject to: ; 1 j N p ; 0 k N ð17cÞ where A, B, C, and D denote matrices as in (10a, b).

Simulation results
Simulations were carried out by using five realizations of 3D RC networks with N X = 5 rows, N Y = 5 columns, N Z = 3 layers, R S = 0.1, R = 1 and C = 0.5 (normalized units). These realizations were obtained by randomly varying the allocation of the R, C components within the network, as in ). The fraction of capacitors with respect to the total number of components (capacitors and resistors) was set to 0.5. From one of the realizations, a reduced-order model of order three was obtained. The model reduction procedure consists of removing dynamic modes which have weak influence on the input-output behaviour of the network. This model, (henceforth termed design model) was used to describe the multi-path propagation process and was the one employed for the calculation of the voltage input sequence. The optimal sequence was then applied to five models associated to different network realizations as a means to verify the robustness of the proposed solution. The imposed constraints (in normalized units) were: |u(k)| B 1.5, |i(k)| B 0.5 and |u(k)i(k)| B 0.125, each of them in normalized units. The sample time was set to T s ¼ 0.010, (also in normalized units), and N was set to 20 in the optimization problem (16), (17a-e). Following the work by Kasper (1982) on equalization of multimode optical fibre systems, the desired output signal adopted in the current simulations was a Gaussian pulse. In all figures, ''Simulation model 1'' is the full-order model that was used to generate the reduced-order design model. The package Multi Parametric Toolbox (MPT) (Kvasnica 2009) was used to manipulate the polytopes and IBM-CPLEX was used to solve the resulting MILP problem. Both were run under the Matlab environment. Figure 3 illustrates the dispersion features of the RC network by using a Gaussianshaped pulse as input voltage. As can be seen in Fig. 3a-c (time domain) and Fig. 3d-f (frequency domain), the output current waveform closely follows the input shape if the pulse is sufficiently wide. However, for narrower pulses (i.e. at higher spectral content), the pulse shape becomes distorted, deviating from the ideal Gaussian shape.
The results in Fig. 3 point to the necessity of a pre-equalization scheme to compensate for the network dispersion in order to obtain an output pulse closer to a Gaussian shape. For this purpose, the optimization problem (16), (17a-e) can be solved by using the desired Gaussian shape as the reference waveform. As a result, the input waveform is calculated so that the output follows the reference as close as possible, given the constraints over the voltage input, current output and power. For illustration, the optimization problem was initially solved without any of these constraints. The results are shown in Fig. 4, in which the Gaussian pulse reference with the shortest period shown in Fig. 3 was discretized to cope with the discretetime nature of the optimization problem. It can be seen from Fig. 4a, b, that the output of the design model follows the reference perfectly in the time and frequency domains. As can be seen in Fig. 4a, there is a phase advance of the current with respect to the voltage, which is coherent with the typical response of an RC network. As for the different simulation models with random allocation of the RC components, the main characteristics of the signal are preserved. Figure 4c depicts the results in the current-voltage plane. In this case, the power In order to cope with the source constraints, the optimization was repeated, now imposing constraints on the voltage, current and power. As can be seen in Fig. 5c, the power constraints are now properly enforced. However, the signal at the output does not follow the reference as closely as before, especially at the peak current demand, as seen in Fig. 5a. As a consequence, the output spectrum in Fig. 5b displays an elevation at high frequencies, as compared to Fig. 4b. In practice, such a distortion could be mitigated by transmitting a signal of smaller amplitude (so that the power constraint becomes inactive) at the cost of a worse signal-to-noise ratio.

Prospects for pre-equalization applications in communications
There is considerable scope to further explore the proposed formulation in a communications setting (Kang et al. 1999;Chen et al. 1999;Vanderveen et al. 1998). Line-of-sight links can assume a simple exponentially decaying channel model whereas diffuse links require a ceiling-bounce channel model. In the time domain, the output yðtÞ from the communications channel is given by: yðtÞ ¼ RhðtÞ Ã xðtÞ þ nðtÞ where hðtÞ is the channel impulse response, R is the detector responsivity, nðtÞ is white Gaussian noise due to the lighting in the room and Ã denotes the convolution operator. In indoors communications  where multi-path propagation can be more pronounced, a ceiling-bounce model (Carruthers and Kahn 1997;Barry et al. 1993;Fernando and Balendran 2005) can be used for the channel model: where uðtÞ is the unit step function and a depends on the room size, and the transmitter and receiver position where L is the height of the ceiling, r is the reflectivity of the ceiling above the transmitter and receiver and c is the speed of light. The parameter a ¼ 2L=c is related to the multi-path rms delay spread D rms by noting that a ¼ 12D rms ð13=11Þ À1=2 . The D rms value can be directly associated to the expression in (3) defining a particular network configuration that mimics the channel response. Furthermore, it is also worth noting that optical wireless systems commonly use simple baseband modulation schemes such as on-off keying (OOK), or pulse position modulation (PPM) (Gfeller and Bapst 1979;Barry 1994). Alternatives that are superior from a power efficiency perspective include digital pulse-interval modulation (DPIM) (Ghassemlooy et al. 1998), dual header pulse-interval modulation (DHPIM) (Aldibbiat et al. 2002) and differential pulse-position modulation (DPPM) (Shiu and Kahn 1999). Furthermore, differential amplitude pulse-position modulation (DAPPM) is also a hybrid modulation technique that has attracted interest (Sethakaset and Gulliver 2005) as it incorporates aspects from both pulse-amplitude modulation (PAM) as well as DPPM. The spectral richness of the chosen input function in our simulations imply that such communications protocols can be supported by adopting a pre-equalization approach. Since demodulation of the encoded signals is commonly performed using hard-threshold decision (HTD), maximum-likelihood sequence detection (MLSD) or zero-forcing decision-feedback equalization (ZF-DFE) modalities, it seems that the piece-wise linearization adopted should not have a detrimental effect in the fidelity of the transmitted information. Generally, when intensity demodulation is performed using a direct detection technique, there should be no problem for compensating from channel distortions. In cases where a heterodyning technique would be used, our formulation would need to be augmented to account explicitly for phase.
Ultra-wide-band communications (Taylor 1994;Win and Scholtz 1998;Vaughan and Scott 1999;Saleh and Valenzuela 1987;Maloney et al. 1990) are defined as communications using very short duration pulse signals of sufficiently high rise-time and falltime having a large ratio of 3-dB bandwidth to the signal's center frequency, typically above 0.25. Since future indoor short-distance communication links are likely to use ultra-wide bandwidth radio or infrared transmitters (referred to as short-pulse or impulse radio systems) embedded in the lighting infrastructure of the buildings, it is also appropriate to also consider propagation of such signals. In this case, the channel model has the form: rðtÞ ¼ P n p n ðt À s n Þ where s n is the delay of the nth propagation path, and p n is the received waveform from each path so that sðtÞ represents a weighted sum of time-shifted versions of the waveform pðtÞ. In this formulation, the channel impulse response would be a function of time and azimuth angle: hðt; hÞ ¼ hðtÞhðhÞ and independent descriptions of the multipath time-of arrival hðtÞ and angle-of-arrival hðhÞ expressions would need to be associated with the RC network structure. A greater sensitivity to frequency-selective fading in multi-path propagation will normally be associated in the UWB implementation of the technique.
It is important to further note that the proposed formulation should also have a wide applications domain in free space or in-fibre propagation using coherent sources (Asif et al. 2013). This includes free-space adaptive cohort secure communication where the generated waveform can take advantage of selective frequency fading to minimize eavesdropping. Pre-equalization of soliton pulses launched into optical fibers and dispersion flattening in non-linear fibres have many applications in the optical fibre communications industry and should be possible with the proposed methodology. Re-casting the optimization procedure in the spectral domain would also enable control of wavelength dependent birefringence in optical fibres. The technique also offers new prospects for other more exotic forms of freespace communications such as through spatial hole burning in the atmosphere where controlled filamentation processes would enable the opening of ad hoc communications channels .
Finally, a further advantage of adopting the proposed RC networks to describe multipath frequency selective fading is the fact that they display transfer functions where fractional Laplace transforms can be used to describe frequency selectivity in a parsimonious manner.
The work also opens the possibility for pre-distortion equalization in communications on the basis of identification using fractional order calculus techniques (Jacyntho et al. 2015).

Prospects for pre-equalization applications in photonic networks
In previous works  we have already discussed that the descriptor based formulation and input/output transfer function derivation associated to the adopted RC networks can be used to model attenuation in different dielectric media (e.g. nano-dielectrics, amorphous or porous materials) which are typically characterized using electrical impedance spectroscopy techniques. Although the optimization procedure was implemented assuming voltages and currents, the methodology is generic and can be implemented in a wide range of propagating media across different parts of the electromagnetic spectrum. Such approach makes the proposed propagation models of use well beyond the confines of existing communications modalities and enables us to take into consideration advances in the femtosecond pulse laser community (Wise et al. 2008;Nathan et al. 2005;Keller 2003;Limpert et al. 2003;Xu et al. 2006;Sardesai et al. 1998). In amplified femtosecond systems for example, spectral management at the input port of the amplifier is necessary to ensure that optical components operate within tolerances and will not be damaged by large peak pulse power transients. This is commonly achieved using hyper-Michelson interferometers, rf-excited acousto-optic modulators or liquid crystal masks in 4f optical systems which are programmed to perform a pre-equalization process. Pulse management is currently performed on a trial-error basis using a spectral phase interferometer or frequency resolved optical gating module which monitors the amplifier output after the excitation waveform has passed through the dispersion-compensated optical components of the system with a spectrometer monitoring the spectral input to the amplifier. Essentially this is a black-box approach which monitors inputs and outputs and inferring the processes in the amplifying medium without providing any understanding of the time domain waveform reflected and refracted through the system optical components. As a result amplitude and phase pulse shaping is always ultra-conservative. With the proposed technique, it is possible to adopt a less conservative approach at the laser input port and more importantly perform pulse management in non-linear optical components, this would have applications in quantum optics (photonic crystal sources of correlated photon pairs to be used in polarization entangled photons applications), high harmonic generation (frequency doubling or tripling crystals), parametric amplification at longer wavelengths, THz generation, four-wave mixing and electromagnetically induced transparency experiments. It should also be of relevance to the advancement of ultra-high power microwave and laser technologies e.g. gyrotron research for plasma excitation and target implosion experiments for fusion as performed at the National Ignition Facility (US).
The formulation also paves the way for pre-equalization of optical pulses coupled to photonic crystal fibres (Reeves et al. 2003;Konorov and Zheltikov 2003;Ouzounov et al. 2005). Photonic crystal fibres offer significant promise for ultra-high bandwidth terabyte communications as well as for the generation of atto-second pulses to perform time resolved X-ray studies (Takashi et al. 2004;Drescher et al. 2001;Christov et al. 1997). One of the main reasons for not widely adopting photonic crystal fibres to generate atto-second pulses through non-linear conversion of carrier envelope stabilized femtosecond pulses is that the fibre output does not remain sufficiently stable over a sufficiently long period required to perform several experiments with the generated atto-second pulses. There is a particularly high sensitivity of the non-linear conversion process to laser beam walk-off.
The assumption of slowly varying dynamics in the RC network adopted in this study should enable the development of a time-varying pre-equalization scheme for dynamic compensation leading to a more stable output.
Future simulations related to the specific case of photonic crystal fibres will concentrate on the application of the modal transmission line theory as has been developed for modelling multilayered periodic media as well as to perform modal analysis of arbitrary shape optical waveguides. In order to achieve this one needs to consider the temporal harmonic electromagnetic fields in a dielectric medium as solutions of the source-free Maxwell equations which can be transformed into a system of differential equations with the electric and magnetic field vectors in the x and y dimensions related to capacitive and inductive responses (Karpisz et al. 2015; with expressions analogous to those found in the telegraph equations. Such approach requires the extension of the current threedimensional RC framework to incorporate inductive components, with the descriptor equations associated with that structure and its boundary conditions converted into state space form. This is a topic of current research by our group. Another interesting domain of applications is in polaritonic integrated circuits (Werley et al. 2010) which currently promise significant potential to advance THz spectroscopy. In these systems, the pulses are generated within the crystal and phonons from crystal edge reflections introduce additional frequency selective attenuation to the input waveform. The pre-equalization procedure proposed can ensure spectral richness of the excitation waveform reaching the sample, accounting for phonon band-gap interference, thus enabling persistent sample excitation.
Finally, it is worth noting that although the formulation discussed in this paper was confined to electrical excitation and was explored within a communications setting, it is straightforward to assume a far more general structure for the adopted RC networks with nodes composed of Gibb's free energy expressions related to acoustic, acousto-optic or chemical processes. Such approach would enable the formulation to provide optimal excitation in acousto-optic crystals or even to more complex physico-chemical systems e.g. ferroelectrics. It may, therefore, be concluded that the proposed formulation is of general interest across the Physical Sciences.
Future refinement of the proposed RC network topology will enable the derived models to account for propagation in complex photonic crystal structures (Datta et al. 1993;Sigalas et al. 1994;Benisty et al. 2002;Kafesaki et al. 2002;Foteinopoulou et al. 2007;Sigalas et al. 1999;Reboud et al. 2010;Vasilantonakis et al. 2012;Tasolamprou et al. 2010). Such developments would have applications in studies of dispersion compensation in fibre based communication systems (Knight et al. 1998;Várallyay et al. 2010). Furthermore, with the addition of non-linear components in the network it would be possible to simulate propagation in metamaterials (Koschny et al. 2005) providing an alternative approach to their modelling as presented in recent work (Liaskos et al. 2015). Finally, it is worth noting that experiments have already been proposed to characterize in both space and time the propagation of an ultrashort pulse in a multimode optical fiber where the optical field amplitude and phase at the output of the fiber can be estimated as a function of delay using time-gated spatial heterodyne interferometry (Rokitski and Fainman 2003). It would therefore be particularly useful to use spectral phase interferometry to experimentally verify the validity of the proposed RC modeling approach to emulate multimode amplitude attenuation and phase delay.

Conclusion
This contribution provides a new way of interpreting multi-path propagation in a dielectric medium as well as in communication channels by assuming an analogy to multi-layer RC networks with randomly allocated resistors and capacitors being responsible for the delay and attenuation properties of the medium. The network responses can describe frequencyselective attenuation encountered in path-dependent propagation. A pre-equalization procedure for a Gaussian pulse input that takes into account the capabilities of the transmission source as well as the transmission properties of the medium is developed. The problem is cast within a Mixed Integer Linear Programming (MILP) optimization framework that uses the developed nominal RC network model, with the excitation waveform customized to generate the desired function in a least squares sense at the receiver. Simulations are carried out with different network realizations in order to evaluate the sensitivity of the solution with respect to changes in the transmission medium associated with multi-path propagation. The proposed approach is of relevance where channel identification is ill-conditioned and equalization techniques are difficult to implement (e.g. displaying fractional order dynamics). It is applicable to both established communication modalities across the EM spectrum as well as emergent indoor communications assuming various modulation protocols or UWB schemes. It is anticipated that such formulations will have a much wider applicability to measurement science, the inverse problem and the photonics modelling community.