FPGA implementation of a 32x32 autocorrelator array for analysis of fast image series

With the evolving technology in CMOS integration, new classes of 2D-imaging detectors have recently become available. In particular, single photon avalanche diode (SPAD) arrays allow detection of single photons at high acquisition rates (\geq 100 kfps), which is about two orders of magnitude higher than with currently available cameras. Here we demonstrate the use of a SPAD array for imaging fluorescence correlation spectroscopy (imFCS), a tool to create 2D maps of the dynamics of fluorescent molecules inside living cells. Time-dependent fluorescence fluctuations, due to fluorophores entering and leaving the observed pixels, are evaluated by means of autocorrelation analysis. The multi-{\tau} correlation algorithm is an appropriate choice, as it does not rely on the full data set to be held in memory. Thus, this algorithm can be efficiently implemented in custom logic. We describe a new implementation for massively parallel multi-{\tau} correlation hardware. Our current implementation can calculate 1024 correlation functions at a resolution of 10{\mu}s in real-time and therefore correlate real-time image streams from high speed single photon cameras with thousands of pixels.


I. INTRODUCTION
Fluorescence correlation spectroscopy 1,2 (FCS) is a powerful experimental technique for measuring the dynamics of fluorescently labeled molecules in solution and also inside living cells. It allows one to determine the particle number, the diffusion coefficient, flow speeds and also photophysical and chemical reaction rates (for an overview, see Ref. 3). In FCS the time trace of the fluorescence intensity fluctuations I(t) inside a small observation volume (usually around 10 −15 l = 1µm 3 ) is monitored. The fluctuations originate from particles entering and leaving the focus, or transitions between states having different quantum yields. Faster dynamics of the fluorescing particles also lead to faster fluctuations, which can be quantified by means of a temporal first-order autocorrelation function (ACF): (1) The ACF usually contains features that are spread over several orders of magnitude in time (nanoseconds to seconds).
The standard FCS setup uses a confocal microscope in combination with single-photon sensitive detectors to a) e-mail: jl@dkfz.de; homepage: http://www.dkfz.de/Macromol acquire the fluorescence time trace I(t) from one focal volume. Then the data are fed into a "correlator" (hardware or software component), which estimates the ACF over a certain dynamic range.
In the accompanying paper by Mocsár et al. 4 , a field programmable gate array (FPGA) implementation of such a correlator circuit for use with a confocal microscope and up to two single photon avalanche diode (SPAD) detectors is presented. In recent years, the availability of fast cameras has triggered the development of different imaging modalities for FCS [5][6][7] , which allow spatial mapping of the diffusion coefficient and other dynamic properties by calculating an ACF for each of potentially many pixels. This requires correlation hardware that can process the data stream very fast, ideally in real time, for all pixels of the image sensor. Here we extend the idea of hardware reuse, presented before 4 , for application to imaging FCS.
Our FPGA-based implementation can calculate 1024 autocorrelation functions in parallel and in real time. The dynamic range of the ACFs is τ min . . . τ max = 10 µs . . . 1 s. As a data acquisition device we use the single photon avalanche diode array (SPAD array) Radhard2 8,9 . From our sensor we read a 32 × 32 pixel frame every ∆t frame = 10 µs where each pixel contains 1 bit of information (no photon or at least one photon in the last ∆t frame ). The design presented here could easily be adapted to other image sensors such as customized complementary metal oxide semiconductor (CMOS) or electron multiplying charge-coupled device (EMCCD) cam-eras.

II. MULTI-τ HARDWARE CORRELATORS
A hardware correlator estimates the ACF in Eq. 1 from a finite sequence of intensity measurements with the number of samples T and the integration time τ min for one sample. When discretizing Eq. 1 with this intensity sequence, care has to be taken not to bias the normalization 1/ I with a given set of lag times τ k ∈ N (in units of τ min , so τ = τ k · τ min ). When the full sequence {I n } n=0...T −1 is available after the measurement, Eq. 3 may be evaluated directly for an arbitrary (also logarithmically spaced) set of lags τ k . This gives an unbiased estimation of the ACF ("direct correlation").
To implement our hardware correlator, we use the multi-τ scheme introduced in Ref. 11, which is also illustrated and compared to a linear implementation in FIG. 1 (for a detailed description, please refer to our accompanying paper Ref. 4). The multi-τ scheme uses a set of S "linear" correlator blocks (FIG. 1(b,c)). The input samples I s,n (n is the same index as in Eq. 2) are summed over increasingly long periods ∆n = m s : with I 0,n = I n .
Each of the linear correlators estimates the ACF at P linearly spaced lags where p = 0 . . . P − 1.
In summary, this results in a quasi-logarithmic spacing of estimatesĝ sym,multi-τ (τ s,p ). The advantage of this multi-τ scheme is its simple implementation in hardware and a large dynamic time range with a reasonable number of channels. Its disadvantage is a systematic error introduced by averaging: As shown in Ref. 12 the estimatorĝ sym,multi-τ (τ s,p ) equals the ideal correlation function g(τ s,p · τ min ) (see Eq. 3) convolved with a triangular kernel with width m s : where * denotes the convolution product and Λ(τ, ∆τ ) = ∆τ − |τ | for |τ | < ∆τ and Λ(τ, ∆τ ) = 0 for |τ | ≥ ∆τ , is the triangular shaped kernel.

III. HARDWARE DESIGN
Here we describe how the hardware reuse scheme in the accompanying paper Ref. 4 can be extended to accommodate many more input channels. Instead of the graphical tool used there, here we employed a low-level hardware description language to gain speed and flexibility. This enables us to fine-tune many parameters of the final design, e.g. operational speed, memory usage, logic resource consumption and routing between logic cells.

A. Single-Pixel Correlator
As shown in FIG. 1, a typical correlator is made up from channels, each corresponding to a certain lag time and consisting of a multiplier, an accumulator and a delay element.
The idea of our implementation is to use one single channel circuit to calculate all channels within one multi-τ correlator. This is possible by serial processing of the lag time channels, since every channel's hardware is identical.
The basic arithmetic operation of one channel is to multiply and accumulate (MAC). Therefore we can map its functionality onto a MAC unit which can be found on most FPGA architectures and which is considerably faster than using generic FPGA logic cells. Only about 100 of these can be found on typical FPGAs, precluding any approach with blocks consisting of several lag channels.
As only one circuit is used to process all channels, we use an internal memory block (block random access memory, BRAM) to store their state (i.e., the content of the accumulator and the delayed signal). We implement this circuit by decomposing it into five steps: To increase performance, these steps are executed in an interleaved manner for four channels simultaneously. This "pipelining" scheme is shown in TAB. I. Our five pipeline steps are compatible with the internal pipelining of common MAC units.
From one block to the next the delay time is doubled (m = 2) in the multi-τ scheme, making the input data rate of block s + 1 half that of block s. Hence we need to execute each block only half as often as its predecessor. Thus, a complete multi-τ correlator can be executed in only twice the run-time ∆t lin of a single linear correlator block: A key requirement is that ∆t lin is at most half the integration time τ min of the input signal I n .
As before 4 , a scheduler guarantees that a linear correlator block s is only executed after its predecessor s − 1 has been executed twice. A counter c = 0, 1, . . . is incremented with every execution of any linear correlator block. The scheduler uses the following relations to determine which linear correlator block s has to be executed at a given counter value c (details see appendix): FIG. 2(a) shows the solution of this relation for c values from 0 to 31. In the binary representation of c for a linear correlator block s patterns are evident that can be used to implement the scheduler efficiently. As shown in TAB. II (for S = 8 linear correlators), correlator block s = 0 is executed whenever the last bit of c is 0 b , correlator s = 1 is executed when the last two bits are 11 b and so forth. This scheme uses only simple comparison operations.
Between two consecutive blocks s − 1 and s, adder circuitry is inserted to sum up two subsequent input signal values I s−1,n−1 and I s−1,n . This is done for both the delayed/local as well as the undelayed/global signal, while they are processed in the pipeline.
The linear correlator as implemented here, together with its scheduler and the summation logic, is called correlation processing element (CorrPE). All channel data and intermediate summation results, the so-called pixel context, are stored in a dual-port BRAM which is associated with the CorrPE.

B. Multi-Pixel Correlator
The CorrPE described above is much faster than required to calculate the ACF for a single pixel in the SPAD array. Hence, we can reuse a single CorrPE for multiple pixels by switching between pixel contexts. This "pixel scheduler" uses a double-buffering strategy. While a CorrPE operates on the current context, the previous context is exchanged with the next context to be processed. The pixel contexts are stored in external background memory (SRAM), because they exceed the capacity of the internal memory. In FIG. 3 we illustrate how we reuse one CorrPE for several pixels in comparison to a naïve implementation using one CorrPE per pixel.
In addition to the channel data, the cycle counter c and an accumulator for the local and the global input signals have to be saved. The latter are used for normalization and cross-correlation.
A single CorrPE is used to process an entire column (ACFs of n y = 32 pixels) of our SPAD array. To handle the full n x × n y array of pixels, we instantiate n x = 32 CorrPEs in parallel. An overview of this scheme is shown in FIG. 4. A "data acquisition" circuit communicates with the SPAD array and provides the image data for the correlators. As the image data are streamed out row by row, and each of the n x CorrPEs is only processing data from one specific context (i.e. a specific row), the remaining pixels have to be buffered in 32 FIFOs (first-in first-out memory buffer) localized in external RAM.
In addition our design contains two USB 2.0 interfaces. One is used to send the raw data stream from the SPAD array to the computer, which allows further data processing. We also use these raw data to verify our correlator design. Via the second interface, intermediate and final results from the correlators are transferred to the host computer. The intermediate results allow implementing a live view of the ongoing calculation of the ACFs.

C. ACF Normalization
The intensity values in the denominator of Eq. 3 are typically obtained from monitor channels M τ k , which accumulate the total photon count at a given lag time τ k . In contrast to other implementations, we use only one monitor M 0 (input signal accumulator) per multi-τ correlator. Symmetric normalization (see also Eq. 3) yields the following: After measuring T samples, the number of samples that have propagated through the correlator to a distinct channel (s, p) is T −τ s,p . Under the assumption that I n is a stationary random process and that T > τ s,p , the perchannel monitors M τs,p can be estimated in the following manner: This normalization and model fitting is done on the host computer, since floating-point arithmetic cannot be implemented efficiently on an FPGA. Although the correlation on the FPGA and the normalization on the host computer can easily be done in real time, the fits in most cases cannot. Usually a single curve fit takes tens of milliseconds, thus fitting all 1024 ACFs would amount to ≥ 10 s additional computing time; but this still allows to display fit parameter maps (images) within an acceptably short delay. In current imaging FCS software systems 13 the data are loaded and correlated on the host computer, which for a measurement of typically 10 s takes ≥ 10 min for 1024 pixels at the frame rate of our sensor.
To show that the estimation in Eq. 10 yields good results, we simulated different correlator types in software. The results for a direct estimation of the ACF using Eq. 3 (green), a multi-τ correlator with a monitor channel per lag (blue) and our estimation (magenta) can be seen in FIG. 5, where the data in (a) and (b) were obtained by correlating the input signal I(t) = 1 + sin(2πt/(1.51 · 10 −4 )) for which the exact ACF is known to be g (theoretical) (τ ) = 1+cos(2πτ /(1.51·10 −4 )) (time t and lags τ are unit free). The data in FIG. 5(c) was created by simulating a T sim = 1 s long FCS experiment with one diffusing species 14 . It was computed with our FCS simulation software described in Ref. 15. Further details on the simulation code are shown in the appendix.
For short lags the estimated ACFs resemble the theoretical curves quite well. Multi-τ correlators have an increased absolute error for longer lags, which is due to the averaging described in Eq. 6. This can be seen especially in the case of the sine wave signal. The multi-τ estimates can still be used for FCS experiments, as here the ACFs usually decay to 1 (white noise) for large lag times, and thus the systematic error drops to zero again (for a detailed discussion of this, see e.g. Ref. 12). For τ T sim /10, multi-τ correlators show additional systematic deviations from the theoretical curve and from the direct estimation, because the channels are not averaged over sufficiently many samples to yield reliable results.
Here the multi-τ implementation with multiple monitor channels performs better due to the better estimation of the normalization factor M τs,p .

D. Crosscorrelation (CCF)
Our design can also calculate cross-correlation: between two input signals I (x) (t) and I (y) (t) at the local J (l) and global J (g) inputs of the CorrPE (see FIG. 1 where both signals are tied to I(t) for ACF calculation). This changes the multi-τ estimator Eq. 9 to: Here the normalization uses the monitors M

IV. PERFORMANCE & IMPLEMENTATION DETAILS
The complete design is implemented in two Virtex-II Pro FPGAs (XC2VP40, Xilinx, http://www.xilinx.com/, San Jose, USA), on a LASP development board 16 . One FPGA is used for data acquisition and line reordering, while the other one is used to implement the correlators. The total resource consumption within the second FPGA is around 80%. Using this hardware platform, there are currently P = 8 channels within each of the S = 14 linear correlator blocks.
In Eq. 7 we showed that the complete multi-τ correlator can be executed within twice the time needed for the first linear correlator block, so we devote half of the execution time to the first linear correlator and the rest to the remaining blocks. Therefore a new input sample can be accepted only once every 2∆t lin . Here ∆t lin = 2P + 3 cycles is the time needed to process a new input sample I n in the first linear correlator block. The 3  To estimate the memory consumption of our design, we first look at the pixel context, which consists of 128 words of 64 bits each. TAB. III shows a detailed memory layout. Sixteen of the upper 32 bits of the raw accumulators store the current value of the delay registers J (l) s,p . The lower 32 bits contain the accumulator G τs,p . The monitor channels are 32 bits each. Data handover between consecutive blocks (accumulated local and global signals) is done via 16 bit-wide memory areas, which is sufficient for S ≤ 16. Since in the later linear correlators the accumulated input signals I s,n are multiplied, the sums G τs,p and also the I s,n can get relatively large and may not fit in the 32 bit memory locations available. However, for our FCS application we estimate that the word sizes used are sufficient, since the per-pixel input event rate is limited. For a deeper discussion of this topic, see Ref. 18.
For our 32 × 32 pixel SPAD array the pixel contexts are stored in 32 · 32 · 128 · 64 bits = 512 KBytes of external SRAM. A second SRAM stores the 256 KBytes used for the FIFOs (2048 entries each) that hold the pixel data until they are processed.
The correlator design is implemented in VHDL (very high speed integrated circuit hardware description language). It can be configured via generics (a VHDL feature). Thus we can reuse our design for different needs (e.g. for different sensor sizes and frame rates). We expect our design to show a significant increase in performance when implemented on newer FPGA generations (e.g. Virtex 5 or Virtex 6 from Xilinx).
To ensure the functional correctness of the designed correlator, a simulation of the hardware was tested using random input data. The outcome was compared to the results of a software implementation of the multi-τ cor-relator using the same data set. Both yielded exactly the same results. The comparison was also done using real experimental data from the SPAD array. Again, both results were identical.
To demonstrate the functionality of the whole system, we tested the design using the SPAD array to record an LED connected to a sine wave generator set to 2.5 kHz. FIG. 6 shows the results of this measurement. A fit to the data recovered the chosen frequency.

V. CONCLUSION
In this paper we presented the implementation of an FPGA-based multi-τ correlator design that can calculate 1024 correlation functions in real time at a minimum lag time of 10 µs. To our knowledge this is the largest number of real-time multi-τ correlators implemented so far in a single device. The minimum lag time of 10 µs in our design is considerably longer than that of currently available hardware correlators (e.g. from ALV GmbH, Langen, Germany or correlator.com, Bridgewater, USA and Ref. 19). However, those are limited to at most 32 auto-correlators.
We use our design to correlate the output of a singlephoton avalanche diode array used as image sensor. This combination will allow us to perform imaging fluorescence correlation spectroscopy at sufficiently fast time scales to resolve even the motion of small molecules in solution and living cells.
Our design is flexible, so beside estimating ACFs, temporal cross-correlation functions of different pixels 20,21 or multiple colors 22 (using spectrally resolved detectors) can also be calculated. By definition, Eq. 8 can be rewritten in the following manner, with C s containing the cycles in which the linear correlator s is processed: It remains to show that C p ∩ C q ∀p, q ∈ N 0 , p = q. Obviously C 0 ∩ C s ∀s ∈ N, k, l ∈ N 0 : Without limitations let r < s; ∀r, s ∈ N 0 ; r, s ≥ 2; k, l ∈ N 0 :

Appendix B: FCS Simulation & Fits
To test different types of autocorrelators, we used our FCS simulation program described in Ref. 15. It simulates the 3D trajectories r i (t), t = 0...T sim , i = 1..K of a set of K fluorescing particles performing a random walk in which each 1D step is drawn from a zero-centered Gaussian distribution with width Here D is the given diffusion coefficient and ∆t sim is the simulation time step. We use a Gaussian approximation for the focal illumination and detection volume: where w xy is the lateral width of the profile, γ = w z /w xy is its aspect ratio and w z is its length. The program then estimates the expected number of photons N phot (t) detected during each simulation time step [t, t + ∆t sim ]: Here N 0 is the maximum number of detected photons per fluorophore and time step, while q f and q det are the quantum efficiencies of fluorescence and detection. To account for the counting statistics in the SPADs, the number of photons N phot (t) actually detected during a simulation time step is calculated by drawing a random number from a Poissonian distribution with average N phot (t).
The time series created by the simulation is then fed into the three different autocorrelator software implementations: (a) direct estimation, (b) multi-τ with one monitor channel and (c) multi-τ with multiple monitor channels. The latter two (b, c) simulate the complete structure of a multi-τ correlator. In correlator (b) we use one single monitor channel, which counts the incoming photons only and then uses Eq. 10 for the normalization. In the correlator (c) we use one monitor channel per correlator lag which counts the photons actually processed by each lag.
The simulation parameters for the data in FIG. 5 are summarized in TAB. IV. The average photon count rate in the simulations was I(t) ≈ 157.4 kHz which is -due to the chosen detection efficiency q det = 0.5 -about a factor of 5 higher than what would be expected from a standard confocal FCS setup. This way the signal to noise ratio is high enough to visualize easily artifacts introduced by the correlator implementation. FIG. 7 shows an example of a time trace N phot (t)/∆t sim generated by our simulation.
The resulting ACFs were fitted to a standard 3D FCS model function 3 : with the average particle number N in the effective focal volume V eff = π 3/2 · γ · w 3 xy and the diffusion decay time τ D = w 2 xy /4D. The fits were performed using a Levenberg-Marquardt least squares fitting routine. FIG. 2. Block scheduling algorithm. The counter c in the first row defines the current cycle. The second row shows the block s that is processed in the corresponding cycle. Further below the order of processing is shown: When a block s has been processed twice, the following block, s + 1, can process the sum of the two preceding s-blocks.
FIG. 3. Comparison of a naïve implementation (a: one correlator per pixel) and an optimized implementation (b: reuse CorrPEs for several pixels) of a multi-pixel multi-τ correlator for a 3 × 3 SPAD array.
FIG. 4. System layout and data path -from data acquisition to correlation. One USB interface is used to stream the raw images, the other for streaming (intermediate) results. The first level cache (L1, double buffered) is used to hold the context of the currently processed pixel. FIFOs for row resorting and for context storage use external memory.  0.151 ms)). The simulations on the right (c,d) were created using an FCS simulation. The top graphs (a,c) are estimates of the autocorrelation function using direct correlation from Eq. 3 (green), a multi-τ correlator with one monitor channel per lag time (blue) and our estimated normalization from Eq. 10 (magenta). Graph (a) also shows the theorectical ACF g (theoretical) (τ ) = 1 + cos(2πτ /(0.151 ms)) for the sine signal (light red). The lower graphs (b,d) show the absolute deviation of the estimates from g (theoretical) (τ ) (b) and from a fit to the curves (d). The parameters resulting from the fit in (c) are the same within < 2.5% for all three estimates.
FIG. 6. Distribution of 992 (gray) ACFs taken by our sensor (first 31 columns only), exposed to a 630 nm LED sinemodulated with a frequency of 2.5 kHz. The correlator was running for 1.2 s (131072 samples at τmin = 10 µs). The gray curves with significantly lower amplitude are due to hot pixels of the SPAD array; since these SPADs fire randomly, the corresponding correlation amplitude is decreased. The median, which tends to be less sensitive towards outliers than the mean, is shown in red. The theoretical model g (theoretical) (τ ) = 1 + A · cos(2πf · τ ) is fitted against the median in the interval [10 µs, 1 ms], and is shown until τ = 2 ms. The fit yields a frequency of (2502 ± 4) Hz. The fanning out of the curve at high values of τ is due to unfilled correlator channels. Compare also FIG. 5(a).