Evaluation of a data acquisition system for SPECT (PET)

Evaluation of a Data Acquisition System for SPECT (PET) P. Stenstrom, A. Rillbert, E Habte, C. Bohm, S.A. Larsson’ Dept. of Physics, Stockholm University, BOX 6730 113 85 Stockholm, Sweden ‘Dept. of Hospital Physics, Karolinska Institute/Hospital, 171 76 Stockholm, Sweden

111. PULSE mOCESSING

Abstract A fully digital data acquisition (DAQ) system [1,2] for a cylindrical (NaI(T1)) SPECT camera has been developed and tested. This system digitizes pulses of the scintillation e vents from the closest (7 or 5) PMTs to the interaction point. The digitized pulse waveforms are processed to obtain estimates of amplitude and time of arrival. We describe various processing methods based on both simulated and real data. The results obtained from a set of three neighboring PMTs using 99””rc was 5.2 ns std. dev. in timing resolution, 12.6%energy resolution and the computing time per pulse about 7 ps.

I. INTRODUCTION Pulse shaping and processing has traditionally been implemented with analog components. With new digital hardware it has become feasible to replace more and more of the front end parts of DAQ systems with digital components, leaving at most the preamplifier as analog parts. Digital signal processing methods for pulse estimation have been developed [3,4,5,6] for various other systems and situations. The goal here has been to evaluate the performance of such methods for a SPECT camera using a NaI(T1) and PMT detector.

11. ACQUISITION SYSTEM The SPECT camera being developed consists of four rows of 18 PMTs arranged in a hexagonal pattern around a cylindrical NaI(T1) scintillator. The DAQ system [I] for the entire camera consists of 24 modules, each responsible for three PMT channels. A module uses three 8 bit free running ADCs to digitize the PMT pulses and a digital le vel-trigger to start acquisition of the pulse w aveforms. A dedicated trigger netw ork allows digital trigger pulses to be sent to nearby modules, forcing acquisition of PMT pulses of neighboring channels. On each module, the three waveforms, a 24 bit sample-clock timestamp and “who triggered me” information are read by a local digital signal processor (DSP). The DSP performs a pulse processing algorithm to e xtract amplitude and accurate timing from the sampled waveforms, and handles the protocol for connecting the module to a high speed Fire wire network. At the other end of the network an event builder collects the detector data from all modules and sorts it into whole e vents for position and energy identification. The DAQ system allows for adjustment of a global sampling clock between 5 and 20 MHz, and to set the follo wing parameters for each module: Trigger level Number of pre-trigger and total samples per pulse Trigger generation and distribution odoff 0-7803-5696-9/00/$10.00(c) 2000 IEEE

A. Optimal pulse processing The optimal pulse processing (OPP) studied in [3,5,6] is an iteration of weighted least squares or fits of a refcrence pulse to a sampled pulse. It is shown in [3] that this gi ves the maximum likelihood estimator of the amplitude A and time of arrival z if the noise in the samples is assumed to be a set of zero-mean jointly normal random variables. How to derive the “weighting function” for comparison against analog pulse processing systems is also described. In [6] it is sho wn that a method using general linear combinations of pulse samples to obtain estimates of A and AT with coefficients chosen to minimize the root mean squared errors (RMSE) of A and AT is equivalent to the x2 method above. Also shown equivalent is the frequency domain approach, where the amplitude and timing of the pulse are obtained con volving the sampled pulse with the impulse response function corresponding to the matched filter of the signal and noise frequency spectra. A benefit of the x2 method is its PO wer to recognize deviant events such as pulse pile-up, with a high x2 value indicating a bad fit.

x2

x2

A brief description of the method is given in the tippendix where as explained, iterations are performed evaluating a parameter vector a = [ A AT BJTfrom the matrix equation a = W(t)s. A, AT and B are the amplitude, time origin and baseline estimates, s is a vector with the N samples of the pulse and Wft)is the weight function for the three parameters, a 3 rows by N columns matrix. W(t) is a function of the time of arrival of the pulse with respect to the sampling comb, W(t) must be repositioned for each iteration, or tabulated in time steps spanning the possible positions of the pulse. The OPP method is almost identical to ha ving a FIR filter working on a digital data flo w, the difference being that here the “filter” is aligned to find the optimal position of the acquired signal waveform. This is done only when there is an acquired signal, saving a large number of calculations.

B. Pulse processing DSP implementation To obtbin a f ast implementation, the OPP method w as implemented in DSP assembler with some modifications to simplify the code. A lookup table of W f t )was used and during iterations only the sign of the AT parameter was used to step forward or backward in the table. This eliminates the amplitude calculation except for the last iteration. Using only the sign of AT is also motivated by the error made when linearizing the fit (see the appendix); the possibility of divergence or oscillations during the iterations is limited. At the zero crossing of AT the iteration is stopped and the amplitude and time calculated.

1332

The trigger of the DAQ system in question is a level trigger, so a number of pre-trigger samples are needed to capture the entire pulse. The number of pre-trigger samples needed arise from a pulse which just reaches the trigger le vel, i.e. the number of samples over T,. This number must then be added to the total number of samples needed to cover both low signals with several pre-trigger samples and high signals with only a few.

An alternative method is to drop the T aylor expansion of the appendix and to perform the fit for the amplitude and baseline only. The so-obtained W(t) consists of tw o columns, d f t l and d ( t ) . Observing the outcome of the scalar product d ( 0 . s for different times, t=At.i, and searching for a maximum amplitude, the amplitude and time can be obtained by interpolation using a 2nd degree polynomial. Denoting the points around the maximum with A,.I,A, (maximum) and A,+1,one gets A+:

A#)=--'

- A: 2

+ A; t + A i , t~ +A+:

(1)

[-1,1],

A

Note that more samples may be needed if noise in the system is correlated over longer time scales, ho wever using too many samples raises the probability of pile up.

A

B. Frequency characteristics

where A: = A ~ - A ~=- w , : s - ~ ~ ; - =~ (8w i - w i - 1 ) 8 is the difference in amplitude between tw o consecutive values of and can itself be implemented as a scalar product AiA... The time and value of the maximum are then given by t,,

= -d- + A ? + 1 , A,,

~A:--A:+~

= Ai

A;

+ tmax-

+ A:+, 4

.

(2)

Furthermore, since AAi 2 0 2 AAi+l if Ai is a maximum, the zerocrossing of can be used to find the peak. Note that the d . s still needs to be calculated once to e valuate the term Ai in (2). This is essentially the same as the OPP method, except that the derivative becomes a difference and is obtained after the x2 fit rather than before.

Since electronic systems tend to pick up lo w frequency noise, from e.g. the PO wer line frequency, an investigation of the sensitivity of the OPP method to this was made. Simulation results with a CR- (RC)2 shaped pulse with T FWHM = 500 ns ) 12 samples sampled at 7 MHz (3.5 samples over T ~ H Mand per pulse are shown in Figure 1. The plot shows the RMSE of

-

- amplitude RMSE ["I/'

time ongin RMSE [TFWHM]

/0.06

IV. NUMERICALSTUDIES A. Sampling frequency and number of samples The first question to address is what sampling frequency is to be used and how many samples per pulse are needed. In the present case the analog pulse from the PMT and amplifiedinedriver can be approximated by an CR-RC shape with FWHM of about 500 ns, b ut in general se vera1 other pulse shapes [7, appendix A] are used in ionizing radiation detectors.

Io.00

1

10

100

1000 frequency [kH2

the amplitude and delay estimate from the OPP method as a The sampling frequency must be at least twice the highest function of the frequency of an added sine signal. The sine signal had a random phase and an amplitude of 5% of the signal. frequency component of the pulse to fulfill the Nyqliist criterion. A common practice is to use f -20dB, the frequency above At each frequency the RMSE of the estimators was calculated which all components are reduced by 20 dB as the highest fre- using 400 simulation runs. The estimators of the OPP method quency. An alternative, fgg%&he frequency below which 99% are unaffected by noise up to about 100 kHz, which means that of the energy of the pulse is contained w as also tested and effectively all mechanical vibration and PO wer line frequency found to give results within 20% of f -20dB Table 1 contains pickup are suppressed. some of the pulses defined in [7] (normalized to their T ~ H M V. EXPERIMENTAL STUDIES and pulse height) with the number of samples needed across TFWHM, To,,, and T, to fulfill the Nyquist criterion. T 0.01 is To evaluate the performance of a DAQ module, several tests defined as the time the pulse spends above 1% of its peak value were done. For these tests the high voltages of the PMTs in the (where dV/dt=O) and T , the peaking time, is the time from system were individually adjusted to get the photo peak at where the pulse passes 1% of its peak value until the peak. about 80% of the ADC range when a collimated gamma source was placed just in front of the PMT. The baseline was adjusted Table 1 Number of samples needed across the pulse for various pulse shapes. to about 2%of the ADC range. This is done to all0 w for drifts in the system affecting amplification and baseline. Samples Samples Samples Pulse Assuming that the worst case TwHM for all PMT channels over To,ol over Tp overT, shape in the full system lies between 400 ns and 600 ns and that the TFWHM To.01 1.29 2.07 shapes of the pulses are CR-(RC)2, the minimum sampling fre2.58 4.15 Gaussian 1.61 quency needed becomes 2.06/400 ns=5.2 MHz. A value of 7.5 MHz was chosen. This means 4.5 samples/rwHM for the 600 CR-(RC) ns pulse and the minimum number of samples become r0.57*4.51+ r2.86*4.51=3+13 samples. The dead time of the 0-7803-5696-9/00/$10.00 ( c ) 2000 IEEE

1333

trigger in the system is 5 samples, limiting the minimum time between two pulses to 2.8 ps or-357 kpulses/second maximum throughput rate.Considering the counting statistics to be Poisson distributed, and-accepting a pile up (events closer than 2.8 1 s ) probability of 5%, a maximum a verage localrate of -In( 1 - 0.05)/ 2.8 ps = 18 kpulses/second. The data presented below was obtained frbm off-line analysis of the pulse waveforms, using &e OPP method and floating point preci-sion in a high-le vel language.-Eor the results presented the full-noise autocorrelation function w as used, although almost no difference was seen when disregarding it.

%

I

[LSB] 2.0 I

I

0.5

Figure 3: Typical Autocorrelation.

A. Acquisition of the reference pulse f and noise covariance matrix V 1 * A collimated source of 99mTcwas placed at the centre in front of the PMT corresponding to the channel to be calibrated. The pulse waveforms acquired were then scanned to find a representative starting reference pulse, f , and the first pre-trigger sample was used as a coarse baseline estimate. Disre garding the noise autocorrelation at this time W can be calculated. Using the OPP method with the W just mentioned on 1000 pulses (rejecting the 5% with highest x2 value as potential pile up pulses) the pulses are resampled to make them appear at the same time as f . Adding them on top of each other and normalizing the result gives a ne w estimate of f. This pulse is then used as the starting reference pulse, and the above procedure is repeated once more. A typical pulse obtained is sho wn in Figure 2.

sample. Note that the v ariance of the time estimate z is inversely proportional to the square of the amplitude.

C. Tests with 241Am, 99mTcand a test pulse. For these tests the different sources were placed in front of the symmetry point between the PMTs belonging one module, and the high voltages adjusted as before. The preamplifiers of the PMTs have an input to allow the application of a test pulse at the same place as the PMT signal. The test pulses were applied using a pulse generator connected to all three inputs.. I ) Energy resolution The energy resolutions (the FWHM of the peak divided by the position of the peak) obtained in the four tests are gi ven in Table 3, and a typical spectrum is shown in Figure 4. counts

0.4

0.2

2

6

4

10

8

12

14

1000

Sample Figure 2: Normalized reference pulse ffrom a typical channel.

1

1

n

12oq

'0

El-\ 200

Using the residues after the noise-free fit ab0 ve gives the remaining effective noise (from the ADC quantization noise, electronic noise and a possible count-rate dependent pile up noise as described in [6]) present in the system. The noise autocorrelation can then be estimated as

0

q 1

e,eJ+,

with i = O...N-1,

(3)

-

,=I

source

where eJ is the residue of the j:tli sample in the fit. The elements of v" are then obtained as V-l,l= Rss(i-j).The autocorrelation of a typical channel i s shown in Figure 3. .

L41Am 'YYC

Test pulse

B. Variance of the estimates With fand V in hand, the estimated parameter CO variance matrix (FCVF)-' (see the appendix) can be formed and is shown in Table 2. aseline B covariancesare expressed in LSBs bf the ADC and the amplitude-time AT in LSBxtime per 0-7803-5696-9/00/$10.00 (c) 2000 IEEE

200

Eneri;&rbitrary scale)

- .

N-I

Rss(i)=

100

50

0

I 1

I 1

Table 3 Energy resolution [%I. channel 1 ]-channel 2 channel 3 *I summed

26.8 21.4 2.7

I

I

30.9 22.2 2.1

I

26 20.3 26

I

17.3 12.4 1.2

2 the moment the true coincidence detect can not be assessed, since detector channels aremot available at two opposite sides within the camera. However the time accuracy from the pulse processing algoiithms can be obtained with

1334

-

the following argument: The estimated time of arrival for chan. ne1 L can be wfitten ~

.I"

c_

~

t, = t+n,

where t is the true time of the e vent and n, a rando withmean p, and-v ariance ssumigg that tde n, are independent, the time difference tG= t, - t j is-formed with mean and variance pt,, = p, = p, - p, and

0:

=

0:

=

0:

+ 0: with dII = 0.

(5)

Only d',,and plr can be obtained from measured data since the true time t is not known. However 02,can be deduced if 3 or more channels are considered, e.g. for 3 channels

acquired pulse wavqform for different number, of samples per pulse. For each curve a minimum can bg: se position of these minimas as 0: = ( G ~ ~ + ( J ~ ~ etc. - G ~ ~ ) / ~ (6) ples per pulse is plotted in Fi be interpreted with care, since lacliing a'model-for the behavior A general time estimate can then be formed as a weighted sum of the curves in Figure 5 a simple 2nd degree polynomial was used to create the figure. of the time estimates of the individual channels involved,

(7) where A, is the weight given to channel i. The mean value and variance of (7) are pT = t

+ pi and

5 4

2 0;

= (CA:O:)/(CA,)'.

(8)

A natural way of choosing A, is to use the estimated amplitude of the corresponding channel, since the v ariance of the time estimate is (section V. B.) inversely proportional to the square of the amplitude. To evaluate the timing accuracy for the tests, a time estimate T3 was formed with equal weight for each of the three channels giving a std. dev. of o;,

6

= (c& + o:, + & ) / 1 8 .

1

0

samples per pulse

As seen from Figure 5 more samples per pulse only provide a better estimate if an appropriate time span is chosen:The increase of the optimal time ;pannea for low a number of sam-

(9)

The time difference std. dev. calculated from time estimates obtained with the OPP method and the resulting std. dev. of the T3 estimate are shown in Table 4. Data is only from events with energy within the photo peak. Table 4 Time difference and T3time estimate standard deviations (ns).

risk for pile up increases and that the signal and noise are already sufficiently acquired. Above about 25 samples the curve drops and no e xplanation apart from the uncertajnty in the data underlying the plot has been found.

E. ADC resolution A further use of the T3estimator is to e valuate how t h l

D.Sampling frequency and numbe Since the energy resolution of a scintil tem is quite high, it is hard to see an y effe (energy) spectrum due to the v ariation of and number of samples per pulse. Ho wever the T3estimator from above-can be used to assessdhe quality of the x2 fit. Figure5 shows the std: dev . of ,T3plotted against the time spanned (samples per p u k e / s-ampling.frequenc' y,) of the

method using the fourier transform for doing this is described in [ 5 ] .

F: Event handling capacity of the system

The DSP implementation of the OPP method finds the zero crossing almost always within 3 to 4 iterations. It tak es the To see the beha vior of the parameters as a function of a DSP about 200 to 300 cycles per pulse to read out, process and pulse shifted in time, the response function is defined as d t ) = format the result into a pack et for transport to the netw ork. W(t,,yAf(tto-t)+BJ=W(tdlAflto)-Atf’(t,,J+Bl for small t. From the Considering that an event has pulses to be processed in all three relation WF = (E)VF)-’FVF = IM the behavior of r(t) around t=O channels of a module, less than in average 20 ps of computing can be examined. r(t)=r(O)+tr’(O) = [A 0 Bi+ + t[O -A 0 7 . So the time needs to be spent per event. amplitude and baseline estimates are at a maximum (or minimum), insensitive to small time shifts, whereas theA.t estimator gives the time origin itself, scaled by -A. VI. APPENDIX Assume that the pulse shape is known apart from its amplitude A, time origin z and baseline B, its sampled expression is g, = g(t0, i, At, A, 7,B ) = Af(to + i A t - 2)+ B ,

(10)

VII. ACKNOWLEDGMENT This work has been supported by a grant from the Knut and Alice Wallenberg Foundation, Stockholm, Sweden.

where fit) is the function describing the reference pulse,

VIII. REFERENCES

ti = to+ idt the sample points, to the start time and dt the time

per sample. T o linearize the dependence on 2, a T aylor expansion around each ti is made; flti-d -fit$Tf’(ti, with f ’ = dfldt, and the function to fit becomes g i = Afi - A z f i + B . This can be represented in a more compact style using column vectors (lower case) and matrices (upper case), e.g. t g = [go ... gN-l] , where stands for transpose. Rewriting the linearized version of (10) using this vector notation gives

[ 11

P. Stenstrom et al. “A N e w Scalable Modular Data Acquisition System for SPECT (PET)” IEEE T rans. Nucl. Sci., Jun 98, 1117-1121.

121

A. Rillbert et al. “A General Firewire Data Acquisition platform Applied to SPECT”, IEEE NSS conf. record 1999.

(1 1)

[3]

where F is a N by M=3 (M