SPAD-Array Contention Signal and Noise Model Suitable for Multilevel Modulation Schemes With Signal Processing

The first signal model for a single photon avalanche diode (SPAD)-array communication receiver for multilevel modulation schemes is reported. This paper proposes a novel, generalised SPAD array signal and noise model for both digital and analogue, synchronous and asynchronous SPAD readout arrays, which includes the competition between the input photons, dark counts and after-pulsing counts. With this contention signal and noise model, multilevel signals including the signal variation after distortion or equalisation can be evaluated. Also, we report the first numerical investigation for SPAD-based, high data rate, free space, visible light communication using higher order pulse amplitude modulation (PAM) with matched filter, linear and non-linear Volterra post-equalization. Simulations have been carried out to analyze and compare the bit error rate (BER) performances under a variety of conditions. The model is verified by comparison with published experimental results.


I. INTRODUCTION
With the ever-increasing demand for data traffic, RF systems are predicted to be unable to support the growing demand for wireless communications. Visible light communication (VLC) is being investigated as a technology to complement WiFi in future heterogeneous communications networks. Currently, for non-optically amplified optical communication links that use intensity modulation with direct detection, the highest optical sensitivity is achieved using avalanche photodiodes (APD). However, the excess noise generated within the APD limits the receiver sensitivity. Operating the APD above its breakdown voltage, as a single photon avalanche diode (SPAD), can eliminate this excess noise. Therefore, SPADs are gaining a growing interest for use in VLC as the most sensitive possible receivers. With the application of technologies such as massive multi-input multi-output transmission, millimeter-wave The associate editor coordinating the review of this manuscript and approving it for publication was San-Liang Lee .
(mm-wave) communication and non-orthogonal multiple access scheme, 5G mobile communication has significantly increased the network capacity [1]. However, the 5G network is still ground based [2]. The 6G network is expected to provide global wireless connectivity from space to underwater. Therefore, SPADs receivers with extremely high sensitivity are potential candidates to complement the terrestrial communication and expand the wireless coverage [3]. Recently, a number of studies using SPADs as VLC receivers have been reported [4]- [12].
A drawback of SPADs is that they need a finite time, typically a few nanoseconds, to recover from each detected photon event. This drawback, associated with a recovery or dead time, can be overcome by implementing an array of SPADs. In addition to the dead time, SPAD arrays have other non-linear characteristics associated with dark counts, afterpulsing, integration period, symbol period, received optical irradiance, array size, photon detection efficiency, data rate and noise [13]. Despite these non-linear characteristics some promising transmission results have been obtained. VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ For instance, recently 500 Mb/s, −46.1dBm sensitivity transmission has been achieved using a 2.8 mm by 2.6 mm custom SPAD array consisting of 4096 SPADs and associated circuitry [10]. With a received irradiance of 1.4 mWm −2 , this receiver operated at 400 Mb/s with a BER of 1.8×10 −3 using on-off-keying (OOK) modulation and decision feedback equalization (DFE). More recently, a commercially available 3.07 mm by 3.07 mm 5676 SPADs array has been incorporated into a VLC receiver. Without any equalization and 500 lux ambient light, this receiver achieved up to 400 Mb/s data rates using OOK modulation [11] and a BER of 10 −3 was obtained with a received irradiance of 0.48 mWm −2 . This is equivalent to a sensitivity of −53.4 dBm which is only 8.7 dB above the Poisson limit. Both receivers are more sensitive than receivers using standard APDs, as reported in a recent study [12]. However, no multilevel BER model suitable for implementing equalisation was reported. Multilevel modulation schemes are normally employed to improve the spectral efficiency of optical communication links to achieve higher data rates [14]. None of the mentioned studies used spectrally-efficient modulation schemes to make the most of the limited SPAD-based link bandwidth. Therefore, in this paper, the first SPAD-array contention signal and noise model suitable for multilevel modulation schemes with signal processing is proposed. The model is verified by comparison with published experimental results. For the first time, using this model, the signal variation after equalisation or high speed distortion can accurately be described. This model thus provides a foundation for simulating future high speed SPAD-based communications.
It is known that bandwidth efficiency can be increased by n times if one symbol is transmitting n binary bits. In that sense, one symbol could have 2 n possible levels. Therefore it is feasible to increase the transmission capacity for the SPAD-based links using higher order PAM modulation. However, the Nyquist condition does not always hold for multilevel modulation unless strong equalisation is used at the receiver. Use of nonlinear equalisation with higher order PAM modulation to increase the transmission capacity of the SPAD link is reported in this work. For the first time, this paper uses the 'contention model' [15] to accurately calculate the signal-to-noise ratio (SNR) of multilevel equalised SPAD links. Moreover, for the first time, BER results of multilevel equalised links are theoretically evaluated and verified by comparison with published experimental results. Last but not least, by simulation, higher order PAM modulation along with matched filter, linear and second order non-linear Volterra post-equalizations are used to show the potential for improving the SPAD link transmission capacity towards gigabit per second data rates.
In Section II, the previously published detector pool model, is reviewed. This is followed in Section III by a description of the new contention model. Section IV lists the contention model novelty. To verify the contention model, simulated signal waveforms using a µLED are compared with published experimental results in Section V. In Section VI, the system simulation model, for a free space SPAD system, is introduced. This is followed, Section VII, by the SNR analysis. Section VIII describes the results of 1 Gbit/s signal equalisation and recovery. BER calculation for higher order PAM SPAD array is derived in Section IX. The novel noise model for SPAD systems is introduced in Section X. Section XI demonstrates nonlinear distortion in SPAD arrays. To verify the proposed signal and noise model, simulation results are compared with published experimental results in Section XII. Section XIII gives the outlook of future SPAD arrays for VLC communication using µLEDs. Conclusions are made in Section XIV.

II. DETECTOR POOL MODEL
This section reviews the detector pool model [16]. Two important models were published by the same research group [17], [18]. In addition, some other models provided further insight in SPAD simulation [19]- [23]. The detector dynamic response can be modelled as a step response to assess the ideal parameters for SPAD arrays. Initially, all the SPADs are available in the array. Once a number of photons are detected, the diodes detecting the photons are under breakdown. A dead time is required for the diodes to recharge back into the pool of SPADs. Therefore, the effective detection efficiency for the following photons of the SPAD array, ePDE, reduces within the recharge time of the SPAD and recovers afterwards. The reduction in detection probability leads to decreased output counts despite continuous input optical power levels.
The detector pool model uses OOK data only and the inter-symbol interference (ISI) is simply due to the step response from the 0 −→ 1 −→ 1 data transition. Theoretically, the model cannot model the ISI if the dead time is greater than the symbol period. Also, the ISI from the 1 −→ 0 −→ 0 and 0 −→ 0 −→ 0 data transitions cannot be modelled. The detector pool model is based on the assumption that the input counts prior to the step change are exactly zero, which is possible in non-return-to-zero (NRZ) modulation but not for higher order modulation schemes. Only first order effects were considered.
For the detector pool model, received output counts C out (n) at a given discrete time n, can be expressed as Equation (2) [16].
where C in (n) is the number of input photons, P is the input optical power (Js −1 ), h is Planck's constant, v is the optical frequency, T I is the integration period (s), ePDE(n) is the effective photon detection efficiency which satisfies the following condition The ePDE 0 in Equation (3) is the initial photon detection efficiency, which is equivalent to the sub pixel photon detection efficiency (PDE). A sub pixel is a single SPAD with its associated circuitry. The instantaneous number of devices that are available, N avail , satisfies the condition N avail = N array , when n = 0. The step response model is a recursive model where the current available detector number N avail is dependent on the number of available detectors in the previous steps. The recursive equation for the step response is given by Equation (4), which is for the case τ d ≤ T I = T , where τ d is the dead time, T I is the integration period and T is the symbol period.
For the detector pool model, the symbol period T was assumed to be exactly the same as the integration period T I . Fig. 1 illustrates a modelled set of step transients, which oscillate about the normalized one level, with increasing optical flux. The whole range from 0 to 1 is not shown in the step response as only the portion near '1' is of interest for signal processing. 'InCounts' is the number of input photons per integration period or per symbol period. In this ideal case, the dead time τ d is assumed to be the same as integration period T I and symbol period T , which here is equal to 1/f RO = 10 ns, where f RO = 100 MHz is the readout frequency. As the detector pool model is a discrete model, the Time axis used in all the figures measures how much time has elapsed after the first integration period, so the first integration is at Time = 0 s. This first order discrete time model fits well with the experimental results, and it replicates the initial peak, the dead time trough and the long-term steady state with good accuracy [16]. The initial peaks are caused by the limited counting resources of SPAD arrays and the associated variation in N avail (n) per Equation (4). It can be seen that increasing the optical flux leads to a stronger nonlinearity. The step response contributes to ISI directly, modifying the receiver total detection efficiency.

III. CONTENTION MODEL STEP RESPONSE
The detector pool model only works in the condition when the integration period is exactly equal to the symbol period and the modulation format is OOK. This section proposes an extension to the model, called the 'contention model', which to first order can simulate the case that the symbol period is not equal to the integration period. In addition, this novel model can model higher order multilevel modulation schemes and can describe the 1 −→ 0 −→ 0 and 0 −→ 0 −→ 0 data transitions. This generalised SPAD array signal model includes the noise from dark counts, after-pulsing and ambient light. Moreover, it is applicable to digital, analogue, synchronous and asynchronous SPAD readout arrays, including the competition between the input photons, dark counts and after-pulsing counts. Lastly, the signal variation after equalization or high speed distortion can be described as well.
A SPAD array can be read out either in synchronous or asynchronous mode. Synchronous mode is normally used to improve the ISI and bit error rate (BER) performance. In synchronous mode, the integration period and the time at which the receiver is optically active, could be used to screen transmitter transitions or optimally place the integration period within the symbol [16]. For asynchronous mode, the integration periods are not synchronised to the signal symbols and BER is normally worse. As a SPAD receiver uses a clock to perform summation, integration, readout and resetting operations, synchronisation between the transmitter and the receiver is critical for good performance. In this work, while in synchronous mode, we assume the SPAD system is perfectly synchronised where a symbol period is a multiple integer of integration periods T = i × T I (i is integer and i ≥ 1).
For a linear time-invariant (LTI) system we normally use the impulse response with unit area to characterise the signal channel without scaling the signal. As most SPAD arrays are non-linear and read out digitally, the step response rather than impulse response is of more interest.

A. SYNCHRONOUS AND ASYNCHRONOUS READOUT STEP RESPONSE
Assuming the transmission format is intensity modulation using ideal pulse amplitude modulation with M levels (PAM-M) and non-return-to-zero (NRZ) pulse shaping, the normalised, continuous time, transmitted signal can be expressed as: where, K is the total number of symbols in the transmitted sequence, T is the symbol period (per the normal notation for PAM), q k is the optical power of the k th transmitted symbol and P(t) is the ideal NRZ pulse response having a value of one for a duration T . This type of signal is a valid model for a µLED transmitter when the symbol period is large compared with the rise-time of the transmitter as is illustrated in Fig. 2a.  If the impulse response of a µLED is exponential denoted as h(t), then the continuous time received signal at the SPAD array input, A(t), can be expressed as: An equivalent idealised received signal, A in (t), can be defined as where, a k is the average power of the k th symbol input to the array During each symbol A in (t) has a constant level. An illustration of A in (t) is plotted on Fig. 2. For a synchronous readout circuit, the readout frequency of the array is denoted as f RO and the integration period T I = 1/f RO is assumed to be the smallest time unit t, i.e. T I = t. The signal symbol rate is SymR and signal symbol period is T . The k th signal symbol is ranging from the n th to (n + i) th integration periods. Based on the pool model in [16], for a SPAD array, the received signal power at discrete time, n, can be expressed as: where a k is the average input signal power of the k th symbol. As Equation (9) is used to calculate the step response, the average input signal power a k for the whole k th symbol period is included in the equation, rather than the average signal level of each integration period, as demonstrated in Fig. 2.
A Out (n) is the average output signal power for the n th integration period within the k th signal symbol period. C a (n) is the average received number of photons per normalised signal power per integration period which when multiplied by a k equals the average received number of photons at time n. C b (n) is the average number of background ambient photons during one integration period.
DCR(n) is the dark counts per second for the whole array. DCR(n)/f RO is the dark counts per integration period and DCR 0 is the initial dark counts per second for the whole array before competition. If the symbol a k were continuously received, then the resulting steady state output counts per integration period is denoted as s k . P ap , is a dimensionless parameter ( 1) which accounts for after pulsing. FF is the SPAD array fill factor (the ratio of photo-sensitive area to total imaging or pixel area). ePDP(n) is the effective photon detection probability of the whole array which measures the probability that a photon, or a dark count, or an after-pulsing count, can be successfully converted to a output count. The effective photon detection efficiency ePDE(n) can be expressed as ePDE(n) = ePDP(n)×FF. ePDE 0 is the initial ePDE which is physically equal to the sub-pixel PDE.
The number of input photons at time n C in (n) and the number of output photons at time n C out (n) can thus be expressed as: so that N avail (n) is a recursive model that depends on the previous available SPADs N avail (n − 1). Assuming the input photon rate is below the maximum average count of the SPAD array and negligible spatial and temporal coincidence, there are 3 different cases to be discussed separately: 1) If τ d < T I , assuming τ d ≈ T I , T I is long enough for the used SPADs to be charged back into the pool, then N avail (n) can be expressed by Equation (4).
can be rewritten as Equation (14) as the SPADs do recharge back into the pool but (m − 1) × T I later than the previous case where m = 1.
3) If τ d = T I , Equation (4) can be rewritten as Equation (15), which represents the ideal case omitting quantization error from readout circuits. This equation represents the ideal asynchronous mode when T can be smaller than τ d to model high speed scenarios. This is a theoretical case where SymR is significant however the large τ d would cause more quantization error.
In the discrete model (Equation (15)), SPADs return to the available pool after one dead time. For a dead time an integer multiple of the integration period, as Equation (14), the number of devices recharged back is C out (n − m − 1). In other words, SPADs do recharge back into the pool, but sometime later than the C out (n − 2), previously used. For example, if τ d = 4 × T I , upon a step, the ePDE would decrease for four integration periods before beginning to increase. This increases the depth of the step response trough. Assuming 1024 SPADs and first integration period giving 150 initial counts, the effective efficiency would follow the series: 100%, 85.3%, 72.85%, 62.2%, 67.7%, 70.3%. This is the main source of amplitude non-linearity.

B. CROSSTALK DISCUSSION
SPAD arrays are known to exhibit crosstalk between pixels. This is due to photons emitted from the avalanche region of a SPAD, triggering detection events in nearby SPADs in the array. An experimentally verified crosstalk modeling approach in a multi-SPAD receiver was presented in [24] where crosstalk was the main noise mechanism. There are many other SPAD array applications where the crosstalk may be neglected where the probability of the crosstalk is often at the 1% to 2% level in SPAD arrays [25]- [28]. From the signal processing perspective, the signal change is negligible when the crosstalk is around 1% to 2%. It is estimated from [16] that 1% to 2% crosstalk has a negligible effect on the BER. In order to fully address the crosstalk for a specific array, a crosstalk threshold needs to be defined in accordance with other noise mechanism levels and the SPAD array performance.

IV. CONTENTION MODEL NOVELTY
The detector pool model can accurately simulate isolated OOK symbols in discrete time, synchronous readout mode. It does not consider signal data rate (not equal to readout frequency), noise, counts contention, ISI (for 1 −→ 0 −→ 0 and 0 −→ 0 −→ 0 data transitions), signal step-response non-scaling for higher order modulation, asynchronous mode and average received optical power. It cannot model ISI and channel dispersion which are pronounced in high speed and long distance transmission. The following sections state how the new model addresses these issues.

A. AVERAGE RECEIVED OPTICAL POWER
The SPAD array response is a function of the average received optical power. A high average received optical power (high input number of photons per second) will use more SPAD resources and hence introduce more non-linearity to the signal response. The term C a (n), which is defined as the average number of received photons due to the received signal optical power per integration period was introduced in Equation (9). C in of Equation (11) includes both received signal and ambient photons. In this way, C a (n) can quantify the intensity of average received optical power in a uniform way for varied signal levels while maintaining the normalised step response and the same average received optical power. C a (n) can be calculated from the input signal intensity waveform. One step response with a C a (n) sweep is demonstrated in Fig. 3. It can be seen that higher average received optical power introduces more non-linearity to the response. The response shown in this figure is smoother compared with the one in Fig. 1 as it uses 64 points rather than 10. VOLUME 9, 2021

B. DATA RATE SIMULATION
As this is a non-linear system, changing the signal data rate changes the shape of the step response. The notation a k in Equation (9) is a function of the data rate. Fig. 4 shows the results for two data rates with normalisation. It can be seen that increasing the data rate, here from 1 Mb/s to 500 Mb/s, distorts the input signal to the SPAD array. The µLED transmitter used in the simulation had a bandwidth of 150 MHz. Compared with the detector pool model, this contention model does not have to assume that the data rate is equal to the readout frequency. The results in Fig. 4 assumes that each symbol has a length of 64 integration periods. The ePDE is not only updated within symbols but between symbols. Therefore, the responses of the symbols are not identical even when the input signal are almost identical as seen from Fig. 4a. For every graph, the first trough is at 0.25 T due to dead time distortion. From Fig. 4a, at the start of the second symbol, N avail is reduced, as a small number of SPADs were used and taken out from the pool during the last integration period of the previous symbol. C a (n) = 80 is very large for illustration purposes.

C. NOISE SIMULATION
The noise notations DCR(n)/f RO and P ap are included in Equation (9) to account for the dark count and after-pulsing count noise as they are non-negligible intrinsic characteristics of a SPAD array. Fig. 5 shows the results for a dark count rate sweep with normalisation. The increase of dark count rates and after-pulsing increases noise and distorts the signal levels and degrades the step response. It is worth noting that the DCR starts from 2.5 × 10 6 Hz per SPAD, which is intentionally set large. The reason is that the readout frequency is f RO = 100 MHz and 2.5 × 10 6 Hz per SPAD DCR contributes 25 photon counts per integration period. The normal value of DCR = 2.5×10 3 Hz per SPAD contributes negligible photons to the step response in this case.

D. ISI SIMULATION
The contention model can account for the ISI. It is clear from Fig. 4 that the ISI can be described in this model, especially the ISI from 1 −→ 0 −→ 0 data transitions. Two examples illustrating this feature are depicted in Fig. 6. Increasing the total number of SPADs in the array not only improves the available number of SPADs N avail , but also improves the ratio between N avail and N array , which in essence increases the effective photon detection efficiency ePDE(n) in Equation (3). The improvement of the N array also reduces the ripple of the signal response. Most importantly, N avail is updated between symbols so that the information of the previous symbols influences the pool of the current symbol. The yellow curves in Fig. 6 demonstrate the pool updates for 11 symbols. The dead time is equal to 0.25 of each symbol period. A trough is observed at the 0.25 symbol period. Though the second symbol is also a digital '1', a large trough is not observed as the detectors in the previous symbol have been recharged back into the pool. The 9th symbol has less N avail compared with the 3rd and 4th symbol as it has higher signal amplitude due to exponential µLED response, though all the three symbols represent digital '0'. The contention model is flexible and suitable for signal processing, especially with complex signal variation.

E. HIGHER ORDER MODULATION AND NON-SCALABLE RESPONSE
a k in Equation (9) may be used to describe higher order modulation. As high input photon counts tend to distort the signal level and introduce more severe signal troughs at the first dead time (high input counts use more SPADs), s k is included to normalise all the counts within each symbol. Fig. 7 shows the example for the normalised SPAD array step responses with PAM 4 modulation. The after-pulsing parameter and dark count rate are chosen to be very large for illustration purposes. The level 1 signal as depicted in Fig. 7a is fluctuating around 0.3 where 0.3 is chosen for illustration purposes only. The actual steady state value for level 1 is determined not only by the noise, but also the impulse response of the LED, channel, ISI, filter response, ambient light, DCR and equalisation. This steady state value can be calculated from the actual signal.
The higher signal levels introduce more non-linearity as C a (n) is the average received photon counts per signal level per integration period. This unified value C a (n) is employed to match the average received optical power. The average signal levels can be calculated from the actual signal.

F. ASYNCHRONOUS READOUT RESPONSE
Equation (15) can model the asynchronous mode when T is smaller than the τ d = T I . If the integration periods are not exactly within the symbol periods, curve fitting and resampling can be used. Then the process of the previous section can then be employed to analyse the signal.

V. DISCRETE SIGNAL WAVEFORM USING A µLED
In this section, received discrete data waveforms are simulated using a 450 nm, 220 MHz µLED. In the simulation, no filters are included to compare the simulated waveforms with the experimental results from [16]. In the simulation, OOK modulation is applied, with an electrical modulation depth of 2V PP (100% electrical extinction ratio), directly onto the µLED. A pseudorandom binary sequence (PRBS) of length of 2 9 − 1 is used. The light from a single µLED having an emission wavelength of 450nm was focused onto the simulated receiver. The transmission distance is 1 m with absolute dark conditions [16]. The dead time is 10 ns and f RO = 100 MHz. The simulation results are compared with the experimental results to validate the contention model.
The µLED used in the experiment has a maximum −3 dB optical bandwidth of 220 MHz if biased at 20 mA [29], [30]. The µLED was biased at 3 mA as stated from [29] and the µLED parameters are the same as reported in [30]. The bandwidth versus current and the frequency response curves for µLEDs of different peak emission wavelength with injected current of 20 mA were reported in [30]. The −3 dB optical bandwidth was estimated to be 50 MHz with 3 mA bias. Therefore, a −3 dB optical bandwidth of 50 MHz was used in the simulation. In the experiment, an average optical power of 20 µW/cm 2 was received at the 2.4 × 2.1 mm 2 ≈ 5 mm 2 Die size, i.e., 1 µW at the Die. A digital eye diagram for 50 Mb/s NRZ PRBS data with a 100 MHz readout frequency, due to the contention model excluding photon shot noise and LED nonlinear transfer characteristics is shown in Fig. 8a. The eye is open due to the low data rate of 50 Mb/s. The limited µLED bandwidth causes ISI. The noise floor is assumed to be 10% of the total counts, which is shown experimentally to be suitable for preliminary design built-in tolerance to noise [29]. As can be seen from Fig. 8, the data transitions from 1 −→ 0 −→ 0 and 0 −→ 0 −→ 0 can be modelled. As can be seen by comparing Fig. 8a and Fig. 8b, the contention model reproduces many of the key features of the experimental results from [29]. In our simulation, 1.25 µW total power at the SPAD array VOLUME 9, 2021 was assumed, which was set equal to 30 photons per SPAD per integration period. The average signal level is 0.5 for NRZ modulation, so the average received photons per SPAD per integration period per level is 60. The slight discrepancy between the simulated and experimental eyes is likely due to the LED non-linearity response and underestimated noise or ambient light. Another possible reason is the underestimated effective photon detection efficiency in the experiment.
Increasing the data rate towards 100 Mb/s while keeping the same readout frequency of 100 MHz leads to a decrease in the opening of the eye. The degraded response is due to the low bandwidth of the µLED, which causes incomplete transitions. It is worth noting that the upper peaks have decreased along with the reduction of the step response steady state. Fig. 9 indicates the contention model achieves results similar to the experimental ones from [29]. Compared with 50 Mb/s scenario, the BER performance is likely to be worse due to the mid level distributions and the increase of '0' level mean. For clarity, the simulations do not include shot noise. If the actual nonlinear L-I characteristics can be obtained and built into the system simulation model, the simulation results could be further improved.

VI. SYSTEM MODEL
In order to evaluate the performance of free space systems using higher order PAM modulation schemes, each element in the system needs to be carefully modelled. The three essential elements in an optical communication system are the optical source, transmission channel and the SPAD-based array receiver. LEDs are the most popular choice as the optical source due to their low cost. In this section, the 'contention model' is used to model the SPAD-array signal [15]. The LEDs have bandwidth limitations which constrain the overall system performance. FFE + DFE and Volterra equalisations are introduced to compensate these bandwidth limitations. The system model block diagram is shown in Fig. 10 where the µLED and SPAD-based array are modelled with exponential and 'contention model' responses respectively. Free space is assumed to be an ideal line-of-sight (LOS) link with a time delayed delta impulse response. For intensity modulation (IM), the optical power of a source is varied in accordance with the characteristics of the modulating signal. Photon counts per second are physically 'power' (J/s) which represent the PAM signal levels. For PAM-n signals, the photon counts per second standard deviation σ n FIGURE 8. Output digital eye diagram for 50 Mb/s NRZ. The simulation used a dead time τ d of 10 ns rather than 13.5 ns as the original paper [29]. The extinction ratio was assumed to be 11. All other parameters are exactly the same as [29], SymR = 50 Mbaud/s, f RO = 100 MHz, N array = 1024, PDE 0 = 20%, P ap = 0.9%, FF = 2.42%, DCR = 7.27 × 10 3 per SPAD, τ d = 1 × T I = 10 ns, oversampling = 256.
(Counts/s) for each PAM level n can represent the respective noise.

A. FEED FORWARD EQUALISER AND DECISION FEEDBACK EQUALISER
The Nyquist criterion specifies a condition where the original modulated data can be recovered without ISI. However, a transmission channel may break the Nyquist criterion even if the transmitted signal is a Nyquist pulse due to the channel not being flat in the frequency domain. To mitigate the ISI, a filter that is able to compensate for the distortion introduced by the channel is required to reshape the signal back to Nyquist. Equalisation techniques can be used to realise this filter.
In order to compensate the channel dispersion, feed forward equalisation (FFE) and decision feedback equalisation (DFE) are employed to remove the pre-cursor ISI and post-cursor ISI. Noise enhancement is always generated as the high frequency components are amplified by the FFE section as it compensates for the pre-cursor channel dispersion. In this work, the tap values are adapted by minimizing the minimum-mean-square-error (MMSE) metric.

FIGURE 9.
Output digital eye diagram for 100 Mb/s NRZ. The simulation uses dead time τ d as 10 ns rather than 13.5 ns as the original paper [29]. The extinction ratio is assumed to be 11. All other parameters are exactly the same as Fig. 8a.

1) FEED FORWARD EQUALISER
The FFE is assumed to be a finite impulse response (FIR) filter. The incoming signal is sampled at the symbol rate and the delayed samples are aligned for a total length of M . The samples are multiplied by corresponding tap values, which have been pre-adapted to give the desired equaliser response. These weighted samples are then summed to estimate the recovered symbols (see Equation (18)).

VOLUME 9, 2021
Assuming an M -tap FFE, the incoming sample sequence S n can be expressed as

2) DECISION FEEDBACK EQUALISER
The decision feedback equaliser (DFE) can be implemented as an IIR or an FIR filter. The decisions made for the received symbols are weighted and fed back to the summing element of the DFE. The available decisions must come from the past, which means the feedback carries information exclusively from the earlier symbols. Therefore, the DFE can recover the post-cursor ISI.
Assuming an N -tap DFE, the decision sequence D n can be expressed as where d[n] is the decision made for the n th symbol of the sequence and the corresponding tap values for the DFE are then the recovered symbolŝ[n] of s[n] from DFE is expressed asŝ This expression has the same form as Equation (18) (22) If the input of DFE is connected to the output of an FFE, Equation (22) can be expressed aŝ Therefore, the input is s[n] and the output symbol isŝ[n − M + 1]. The combined FFE and DFE schematic is shown in Fig. 11.  Fig. 12 illustrates low complexity linear and nonlinear equalisation employed together. The equalizer consists of a symbol spaced feedforward (FFE) and a decision feedback equalizer (DFE). As depicted in Fig. 12, the solid lines represent the linear equalisation while the solid and dashed lines together illustrate the nonlinear equalisation. The nonlinear equalisation employed here is 2 nd -order Volterra kernels, see Equations (24). Volterra kernels of orders higher than two are ignored due to the dramatically increased complexity which constrains practical implementation [31]. The recovered symbol at time t, x [t], due to the combined linear and non-linear equalization is given by Equation (24): where M and N are the number of FFE and DFE taps respectively, x and x are input and output respectively. The first term in Equation (24)

VII. SNR ANALYSIS
The square root of the dark count rate (DCR) is normally considered as the noise floor [13]. In [32] the dynamic range of a SPAD is defined as Equation (25).
DR dB = 20 log 10 C max √ DCR (25) This assumes that the noise floor is static with increasing flux. Even at the maximum counting rate C max , the dark count can still be detected [32]. In [33], the noise floor σ = √ Photons + DarkCounts + AfterPulses is shown decreasing with optical power as the count rates starts to saturate and deviate from its ideal linear fit. At full saturation, no extra counts above the maximum can be detected. As the photon flux increases, the significant contribution to the noise floor changes from the DCR to the photon shot noise,

√
Photons. SPAD-based array receivers directly convert the received optical signal power (joule per second) digitally to the number of photon counts per second, therefore, the SNR for SPAD-based array can be written as Equation (26). It is worth mentioning that C Signal here is the photon counts purely from the average received optical power, which does not include the noise counts from dark counts, after-pulsing counts, or ambient light counts. In this way, the SNR introduced here can be made consistent with the normal form used for the APD SNR. (26) For Gaussian white noise the BER can be estimated from the SNR assuming the channel is linear and the receiver has sufficient dynamic range and is not saturated nor signal starved. Also, whilst equalisers increase the distance between PAM levels they may also cause some noise enhancement and so SNR can be a good metric for some equalised receivers.

VIII. SIGNAL EQUALISATION AND RECOVERY
A. LOW SPEED 75 MBIT/S TRANSMISSION Low speed data transmission simulations are carried out based on PAM-8 for a back-to-back link. A 2 13 − 1 PRBS sequence is used to generate 1000 PAM-8 symbols. For this scenario, the SPAD signals with and without matched filter, FFE + DFE equalized signals and their corresponding eye diagrams are shown in Fig. 13. From Fig. 13b and Fig. 13d, it can be noticed that the matched filter adds ISI. However, from Fig. 13f, one can see that at the sampling time the level thickness is reduced compared to the raw SPAD signal or the matched filtered signal. Therefore, the SNR is improved with FFE + DFE equalisation. It can be seen from Fig. 13d and  Fig. 13f that linear FFE + DFE equalisation only slightly improves the SNR. This is due to low ISI at low speed. It can also be concluded that the back to back system can achieve 75 Mbit/s transmission without equalisation. FFE + DFE equalisation with a T-spaced 1-tap FFE and a T-spaced 1-tap DFE are used in these calculations.

B. HIGH SPEED 1 GBIT/S TRANSMISSION
This section will take an example to demonstrate the significance of this model for high speed SPAD communication.
Matched filters are employed after the SPAD array response as the optimum receiver filters for the SPAD link. Both linear and nonlinear equalisation techniques are employed after the matched filters to compensate the dispersion. In classic digital communication theory it has been shown that, if the noise is white, then the optimum receiver filter is called a 'matched filter' (MF) [34]. In the absence of noise, if the output of the channel is assumed to have a response to a single isolated symbol equal to h(t), then the optimum receiver filter is a MF having a pulse response equal to h(−t), i.e. the MF has the time reversed pulse response of the isolated pulse response of the channel [34]. The MF receiver is optimum in the sense that it maximises the SNR at the decision instant and it is also an optimum front-end filter for equalizing receivers. In addition, it has been proven that when the optical signal has Poisson statistics (SPAD signal) and the modulation scheme is PAM then the optimum equalizing receiver is a matched filter followed by a transversal filter [35]. Thus, theoretically, equalization could be applied to SPAD-based array links to further improve the performance.
The case of 1 Gbit/s high speed and high received optical power link is discussed in this section. As an initial example, the average received photon counts per signal level per integration period will be set equal to 20 (C a (n) = 20). The SPAD PAM-16 signals with and without matched filter, with Volterra equalisation and their corresponding eye-diagrams are depicted in Fig. 14. Without equalisation, the system cannot achieve gigabit/s transmission even with the backto-back situation due to the severe distortion of the signal introduced by the limited µLED bandwidth and SPADs array intrinsic noise and nonlinearity. Nevertheless, the ISI can be removed if the received signal is reshaped to satisfy the Nyquist criterion by Volterra equalisation. Volterra equalisation with T-spaced 2-tap FFE and T-spaced 2-tap DFE is used in these calculations. The recovered eyes in Fig. 14f are open whereas those of Fig. 14b are closed. This shows the potential for gigabit/s transmission using SPAD-based receivers with higher order PAM modulation and equalisation.

IX. BER CALCULATION FOR PAM-2 m SPAD ARRAY
A simple first order PAM noise model was proposed in [36] where the model assumes the noise in the SPAD array is photon shot noise limited and the noise perturbations that moves between more than one symbol are negligible. The BER model was simplified by the use of Poisson distributions and the BER models for simple PAM-2 were as follows [36].
For a PAM-2 modulation format, the received average photon counts for each of the levels were defined as n 1 and n 2 . At high count rate, the model described in [36] simplifies to a Gaussian distribution with mean and variance both equal to n of the Poisson distribution, the optimal decoding threshold th can be determined analytically by simply finding the intercept between the two distributions [36]: However, Gaussian distribution BER model for SPADbased PAM-n will be assumed since the large number of photons allows the Poisson distribution to approach the Normal distribution. Equation (29) is used to calculate the decision thresholds rather than the half way thresholds between symbols. After equalisation, the received output photon counts can be expressed as: where a indicates the different PAM symbols recovered at the receiver, S a (t) and σ a (t) correspond to the sampled value and noise of the symbol a. If a given PAM symbol is continuously transmitted then the sampled output C out (t) is assumed to be a Gaussian probability density function (PDF) with a mean of S a , and the PDF can be expressed as Equation (31) f where C a is an output signal sample for symbol a. Assuming all the symbols are equally possible, then The total symbol error ratio is SER Total = SER 1 + SER 2 + SER 3 + .. + SER n (33) As illustrated in Fig. 15, the symbol-error-rate (SER) can be divided into three parts, including the first symbol S 1 (leftmost side), the last symbol S n (rightmost side) and the symbols S a (a = 2, 3, . . . , n − 1) in the middle.
For a PAM-(n = 2 m ) system, m is integer and (m ≥ 1, n ≥ a ≥ 1), if a = 1, SER 1 can be expressed as V Th1 is the decision threshold level between S 1 and its adjacent symbol S 2 , and Therefore, Equation (34) can be rearranged to If 2 ≤ a ≤ n − 1, SER a can be expressed as If a = 2 m = n, SER n can be expressed as V Th(n−1) is the decision threshold level between S n−1 and S n assuming the symbols are fixed. As the noise is signal dependent due to SPADs array non-linear response, the decision thresholds are not equidistant from adjacent symbols. Given the threshold expression in Equation (29), the V Th(n−1) can be deduced as follows Assuming natural coding, a simple conversion between SER total and BER total is BER total = SER total m (40)

X. SIGNAL AND NOISE MODEL
The model developed in this section is based on the assumption that the output photon counts have not reached the total counting limit of the SPAD system. All the results obtained are well before the saturation of the SPAD devices. SPAD-based receivers are photon counting detectors, which can convert the optical signal directly into a digital signal as the number of photon counts per second. The average received photon counts N average due to the optical source is expressed as follows where P average is the average received optical power from the optical source, t is the time during which the photons are collected. h is the Planck constant, c is the speed of light in vacuum and λ is the photon's wavelength.
The average received photon counts N ambient due to ambient light is expressed as follows where P ambient is the average received optical power from ambient light. The total average received optical power P can be expressed as Hence, based on the 'contention model' [15], the average output photon counts of the system C out can be calculated as follows where, N DCR is the number of dark count photons during time t. K (P total ), is the ratio of the long run steady state output photon counts over the initial output photon counts after t, which is a function of the total average received optical power P total . t is usually set to one second because C out is usually expressed in photon counts per second. The term K (P total ) can convert the initial output photon counts to the average output photon counts. K (P total ) takes into account not only the actual signal but also the introduced nonlinearity. As can be seen from Fig. 16, K (P total ) decreases as P total increases which describes the SPADs array signal nonlinearity especially in the high received optical power regime.
The average output photon counts C noisefree of the system purely due to the expected average received optical power without noise can be expressed as Only monochromatic light is considered. Assuming the calculated average levels for the PAM signal is L and the photon levels to be equidistant, the average output photon counts per level without noise can be expressed as follows Equation (46) can be used to calculate V Th1 − S 1 , V Tha − S a , −V Th(a−1) + S a , and −V Th(n−1) + S n in Equation (36), Equation (37), and Equation (38) respectively to calculate the BER. If noise is added, the average output photon counts per level with noise can be expressed as For the SPAD array, the noise is photon shot noise limited as the inherent dark count rate noise and after-pulsing count noise in a SPAD receiver is Poisson-distributed due to the nature of the receiver [36]. It would be simplistic to assume that each signal level has the same noise. However, to simplify the calculation, the standard deviation σ of the noise for each signal level can be calculated as follows.
where, as before, a denotes the signal level. C 1 , the counts for the lowest level, is associated with the 'all-zeros' transmission state, should be equal to or above the SPAD dark count rate and ambient light which sets the minimum intensity level [36]. The noise for SPAD sensor can be expressed as where S 1 is the mean signal values of all the symbols received at level 1 with noise. C 1 symbols, which are customarily associated with the 'all-zeros' transmisson state, cannot always be 0 due to the signal distortion through the channel. Knowing the SPAD noise σ n and the average output photon counts per level C perlevel and C perlevelnoise , the SPAD BER can be accurately calculated before and after equalisation. A typical PAM-4 histogram with Poisson distribution approximated by Gaussian with symbol rate of 1 Mbaud is illustrated in Fig. 17.

XI. NONLINEAR DISTORTION IN SPAD
Asymmetrically clipped bit-error ratio performances of SPAD-based arrays were demonstrated in [37]. In the simulation, the maximum optical irradiance and the low error area were defined as the metrics of the nonlinear distortion [37]. In order to generate received PAM symbols, the output of the SPAD array is counted during a time, T I , at time instances t = kT of the received optical signal. T which denotes the symbol period of the time domain PAM signal, is assumed to be longer than the integration period T I . These photon counts are denoted by y(k) which is the sum of the photon counts from each individual SPADs, x i (k). N array is the number of SPAD devices in the array. As the photon counts from each individual SPAD can be estimated using Poisson statistics, the actual photon counts at the output of the SPAD array can be described as where the average photon counts C out (k) (K (P total ) = 1) over the short-time averaged period, T I , at time instances t = kT of the received optical signal was described in Equation (44). The output of the SPAD array is the number of photons, y(k), and the system requires the amplitude of the electrical signal (optical power) to demodulate the received signal to the original encoded bits. Therefore, a photon-to-amplitude equalizer is employed to convert the received photon number to the corresponding electrical signal amplitude (optical power), which can then be scaled to the original signal. In the SPAD-based system, a special form of nonlinear distortion which is caused by the saturation of SPAD devices should be considered. The forms of nonlinear distortion can be classified as passive quenching (PQ) or active quenching (AQ). For the PQ SPADs, any counts including signal, dark count and after pulse occurring during the dead time are not counted, but the dead time is extended. For AQ SPADs, any additional events arrive during the dead time are not registered and do not extend the dead time at the expense of more complex configuration, more area and power. AQ SPADs are non-paralyzable detectors and can potentially achieve higher count rates than PQ SPADs. Table 1 summarizes the parameters used for PQ and AQ SPADs comparison. Fig. 18 shows the average outputs of the PQ SPAD array and the AQ SPAD array as a function of the average received optical power when N array is changed to 4096 and dead time τ d is changed to 1 ns with other parameters unchanged. As increasing the N array or decreasing the τ d , the maximum photon rate improves, and the nonlinear distortion points shift rightwards, extending the linear region. The results using our model are in a good agreement with the results published in [37].

XII. SIMULATION RESULTS COMPARED WITH PUBLISHED EXPERIMENTAL RESULTS
For the following sections, the simulation parameters used are exactly the same as [10], which are listed in Table 2.
The 4096 receiver element array divides into two sets of 64-row XOR trees feeding digital readout chains placed on the left and right of the active area, one is being read out and one is counting. Therefore, the effective number of SPADs is 32 × 64 = 2048. Taking the ambient light and the effective number of SPADs into consideration, the simulated photon transfer curve at 400 MHz RX sampling rate with no dead time nonlinearity, is depicted in Fig. 19a. The results from [10] is also plotted in Fig. 19b for comparison. The simulation achieves a good agreement with the results presented in [10]. The slightly overestimated maximum photon count rate may be due to the neglection of the dead time pile-up distortion and SPAD detector redundancy.  The BER results for 500 Mb/s PAM-4 are illustrated in Fig. 20, the results are capped before the nonlinear distortion region for better illustration. The dead time is assumed to be equal to the symbol period to simplify the calculation. The simulation achieves comparable BER performance with the results published in [10]. Volterra equalisation improves the BER from 2.5 × 10 −3 to 1 × 10 −4 at −46.1 dBm sensitivity. The original experimental results in [10] achieved a −46.1dBm sensitivity at a BER of 7.6 × 10 −3 without equalisation, which is comparable with the results 2.5 × 10 −3 using our signal and noise model. The original equalisation in [10] slightly improved BER from 7.6 × 10 −3 to 2 × 10 −3 .  Our signal and noise model achieves better BER results. Also, Volterra equalisation outperforms the reported equalisation in [10] with only 2 feedforward and 2 feedback taps. In addition, the noise floor is not considered in the simulation which may cause a slightly improved BER performance compared with the experimental ones. The signal model used in this calculation is the one described in [15] with equidistant signal levels. The noise is assumed to be Poisson-distributed for each signal level. In addition, the Poisson threshold (Equation (39)) is employed so that the optimum BER performance is calculated. Moreover, as the ambient light and laser diode collimation details were not reported, the ambient light may be underestimated. In addition, at high photon counts, the dead time changes due to some counts occuring during the dead time (including signal, dark count and after pulse) are not registered, but are extending the dead time. The assumed smaller dead time of our simulation also causes the improved performance. Also, our slightly overestimated photon transfer curve might cause improved BER performance as well. It is worth noting that as the average received optical power increases towards the saturation region of the transfer curve, Volterra equalisation starts playing an important role to reshape the signal back to the Nyquist. Therefore, Volterra equalisation can improve the transmission capacity of the SPADs array especially in the high nonlinearity region.

XIII. FUTURE SPADS ARRAY FOR VLC COMMUNICATON
The results of the previous subsections are based on the SPAD device reported in [10]. The results of this subsection are based on a future SPADs array, which has the parameters listed in Table 3. Assuming ambient light noise and collimation profile is the same as [10], the BER results using the AQ SPADs array are shown in Fig. 21. The Volterra equalisation with 2-tap FFE and 2-tap DFE can achieve BER 1×10 −12 and −25 dBm sensitivity transmission in a practical, background insensitive VLC link at 1 m in 1 klx ambient conditions using a 450 nm µLED. The matched filter slightly improves the BER performance. The increased number of SPADs and reduced dead time significantly prolong the linear region. In addition, the large number of SPADs reduce the signal ripples and improve the steady state value which improves the BER. Volterra equalisation is required when the optical power is high due to the increased second order nonlinearity introduced by the SPAD signal.

XIV. CONCLUSION
A SPAD-array contention signal and noise model is proposed, which is suitable for multilevel modulation schemes with signal processing. The improved model can simulate arbitrary data rates with advanced modulation formats. The model also enables matched filters and equalization to be included in simulations. The proposed model can be used for BER calculations.
To verify this model, the modelled digital eye diagrams for 50 Mb/s and 100 Mb/s NRZ modulation are compared with published experimental results using the same parameters and good agreement was achieved. With this new contention model, the signal variation after equalisation or high speed distortion can be described.
The first numerical investigation for SPAD-based high speed visible light communication systems using higher order PAM with matched filter, linear and second order non-linear Volterra post-equalization is reported. Also, a simple, generalised digital SPAD receiver noise model for higher order PAM modulation has been proposed and employed for the BER calculation of the equalised high speed PAM signals. The use of both Volterra and FFE + DFE equalisations are proposed for SPAD-based data transmission in order to mitigate the non-linear response of SPAD signal. The generalised digital SPAD receiver noise model for higher order PAM modulation and the BER model are also verified by comparison with the published experimental results.