Calculation of time resolution of the J-PET tomograph using the Kernel Density Estimation

TL;DR

This paper presents a novel analytical approach combining Tikhonov regularization and Kernel Density Estimation to accurately calculate the time resolution of the J-PET plastic scintillator scanner, validated through experimental data.

Contribution

It introduces a new method for estimating time resolution using closed-form formulas and validates it with experimental data from the J-PET scanner.

Findings

Analytical formulas for time resolution were derived.

Experimental results agree with theoretical predictions.

Method improves understanding of J-PET timing performance.

Abstract

In this paper we estimate the time resolution of the J-PET scanner built from plastic scintillators. We incorporate the method of signal processing using the Tikhonov regularization framework and the Kernel Density Estimation method. We obtain simple, closed-form analytical formulas for time resolutions. The proposed method is validated using signals registered by means of the single detection unit of the J-PET tomograph built out from 30 cm long plastic scintillator strip. It is shown that the experimental and theoretical results, obtained for the J-PET scanner equipped with vacuum tube photomultipliers, are consistent.

Figures

Tables

Table 1. Table 1: Summary of theoretical CRT calculations of the J-PET scanner. The parameters are described in the text. The presented values of CRT take into account an additional smearing of the time due to the unknown depth of interaction in a scintillator strip with a thickness of 19 mm, see Eq. ( 21 ) for details.

Parameter	Unit	Photomultiplier type
		PMT	MCP
$N_p$	1	220	220	350	700
${\sigma }_p$	ps	300	500	420	360
${\sigma }_{p\left(opt\right)}$	ps	500	500	420	360
Tr$\left(S_p\right)$	V2	0.73	0.55	0.39	0.23
Tr$\left(S_r\right)$	V2	0.22	0.0	0.0	0.0
Tr$\left(S\right)$	V2	0.95	0.55	0.39	0.23
CRT	ps	290	260	215	170

Equations

$$t_r=t_e+t_p+t_d,$$

$$f_{t_r}(t)=(f_{t_e}\ast f_{t_p}\ast f_{t_d})(t),t>0.$$

$$f_{t_e}(t)={\kappa }_e{\int }_{\Theta }^t\left(e^{-\frac{t-\tau }{{\tau }_d}}-e^{-\frac{t-\tau }{{\tau }_r}}\right)e^{-\frac{(\tau -\Theta -2.5{\sigma }_e{)}^2}{2{\sigma }_e^2}}d\tau ,$$

$$t_e⩾\Theta .$$

$$f_{t_p}(t)=\frac{{\kappa }_p\cdot x}{t^2}\cdot e^{-{\mu }_{\text{eff}}\cdot v\cdot t},$$

$$t_p⩾\frac{x}{v},$$

$$f_{t_d}(t)=\frac{1}{\sqrt{2\pi }{\sigma }_d}\exp \left(-\frac{(t-T_d{)}^2}{{\sigma }_d^2}\right).$$

$$t_r⩾\Theta +\frac{x}{v}.$$

$$t^{(n)}=n{T}_s+\Theta +\frac{x}{v}n=1,2,...,N,$$

$$y(n)=\beta (E,x)\cdot f_n,\text{where}{f}_n=f_{t_r}(t^{(n)})n=1,2,...,N,$$

$$\beta (E,x)={\beta }_E\cdot {\beta }_x.$$

$${\hat{y}}_{(L)}({\Theta }^0,x^0)$$

$${\hat{y}}_{(R)}({\Theta }^0,x^0)$$

$$v_{(L)},v_{(R)}\sim \mathcal{N}(0,S),$$

$$\Delta \Theta$$

$$\Delta x$$

$$W(\Delta \Theta ,\Delta x)=(y_{(L)}-{\hat{y}}_{(L)})(y_{(L)}-{\hat{y}}_{(L)}{)}^T+(y_{(R)}-{\hat{y}}_{(R)})(y_{(R)}-{\hat{y}}_{(R)}{)}^T.$$

$$(\Delta \hat{\Theta },\Delta \hat{x})=\arg \min W(\Delta \Theta ,\Delta x)$$

$$\Delta \hat{\Theta }\sim \mathcal{N}(0,{\sigma }_{\Theta }^2),$$

$$W(0,0)$$

$$=\sum _{n=1}^N{v}_{(L)}^2(n)+v_{(R)}^2(n).$$

$$\text{Var}[W_{\text{min}}]=2\sum _{n=1}^{N}2S^2(n,n).$$

$$W(\Delta \hat{\Theta },0)$$

$$\approx {\alpha }_0+{\alpha }_1\Delta \hat{\Theta }+{\alpha }_2\Delta {\hat{\Theta }}^2.$$

$$W(\Delta \hat{\Theta },0)\approx {\alpha }_2\Delta {\hat{\Theta }}^2.$$

$$\text{Var}[W_{\text{min}}]\approx 2{{\alpha }_2}^2{{\sigma }_{\Theta }}^4.$$

$${\sigma }_{\Theta }=\sqrt[4]{\frac{2{\sum }_{n=1}^N{S}^2(n,n)}{{{\alpha }_2}^2}}.$$

$$\text{CRT}=\sqrt{11.04\cdot {\sigma }_{\Theta }^2+\frac{d^2}{c^2}}.$$

$$v=v_p+v_r,$$

$$S=S_p+S_r.$$

$$\tilde{y}=y+v_p.$$

$$\tilde{y}=\sum _{k=1}^{N_p}\widetilde{y_k}.$$

Full text

• 07 April 2017

Keywords: Positron Emission Tomography, Time Resolution, Kernel Density Estimation

1 Introduction

The Jagiellonian PET (J-PET) Collaboration constructs a PET scanner from plastic strips forming the barrel (Moskal et al2011, 2014a). An example of arrangement of the scintillator strips in the J-PET tomograph is visualized in Fig. 1. The proposed setup permits to use more than one detection layer thus increasing the efficiency of $\gamma$ photon registration (Moskal et al2016). A single detection module consists of a long scintillator strip and a pair of photomultipliers attached to the opposite ends of the strip. Measurement with such detector results in timestamps from both sides of each scintillator, which allow to extract the timing, position and energy information of each $\gamma$ photon interaction. The time and position of the $\gamma$ photon interaction in the scintillator strip is calculated based on times at left ($t_{(L)}$) and right ($t_{(R)}$) side of the strip. In the first approximation, the time of interaction may be estimated as an arithmetic mean of $t_{(L)}$ and $t_{(R)}$ and the position of interaction along the strip may be calculated as $(t_{(L)}-t_{(R)})v/2$, where $v$ denotes the speed of light signals in the scintillator strip. The energy deposited in the scintillator strip may be expressed in terms of the number of photoelectrons registered by the photomultipliers and is proportional to the arithmetic mean of a number of photoelectrons registered at the left and right sides of the scintillator; the value of energy calibration factor was evaluated in Ref. (Moskal et al2014b). The registration of single event of positron emission, used for the image reconstruction, is based on the detection of both $\gamma$ photons in two modules in a narrow time window. Therefore, a single image-building event includes information about four times of light signals arrival to the left and right ends of the two modules that register in coincidence. The J-PET detector offers the Time of Flight (TOF) resolution competitive to existing solutions (Humm et al2003, Townsend et al2004, Karp et al2008, Conti 2009, Conti 2011, Słomka et al2016), due to the fast plastic scintillators and dedicated electronics allowing for sampling in the voltage domain of signals with durations of few nanoseconds (Pałka et al2014).

Recently, the time resolution, defined hereafter as the standard deviation, of about 80 ps has been achieved for the registration of $\gamma$ photon in 30 cm long scintillator strips read out at both ends by the vacuum tube photomultipliers (Moskal et al2014b, Raczynski et al2014). Such resolution results in coincidence resolving time (CRT) of about 275 ps as shown in Ref. (Moskal et al2016). Further improvement of time resolution requires developments in techniques of signal processing and effective parametrizations of detector’s features. Our estimate of the time resolution is based on statistical properties of the signals in plastic scintillators. Distribution of the time of the photon emission followed by its interaction in plastic scintillators was described in Refs. (Moszynski, Bengtson 1977, 1979). Following the time order statistics analysis described e.g. in Refs. (Seifert et al2012, Degrot 1986, Spanoudaki, Levin 2011), the statistical framework allowing for the analysis of photon propagation in the scintillator strips was proposed in Ref. (Moskal et al2016).

In this paper we propose a novel approach to calculate the time resolution of the PET scanner based on ideas from the Tikhonov regularization (Tikhonov 1963, 1977) and Kernel Density Estimation (Parzen 1962, Rosenblat 1956) methods. We investigate the quality of estimation of time resolution based on the scheme with a single scintillator strip detector introduced in Refs. (Moskal et al2014b, Raczynski et al2014). The most important aspect of the time resolution evaluation involves the statistical description of noise. The noise in the measured signal comprises two components: statistical fluctuations of the number of photoelectrons registered by the photosensor, and effect of the limited number of samples of the signal in the voltage domain. In Ref. (Raczynski et al2015a), the formula for calculations of the signal recovery error was introduced and proven. In this paper we determine dependence of the signal estimation error on the number and shape of registered photoelectron signals. Theoretical results are compared to the experimental resolutions achievable using traditional readout with the vacuum tube photomultipliers. The method is verified by setting in calculations the same conditions as in the experiment, as described in Refs. (Moskal et al2014b, Raczynski et al2014).

The J-PET tomograph can be equipped with various types of photomultipliers: the vacuum tube photomultipliers (standard in the J-PET prototype), the silicon or the microchannel plates photomultipliers. In case of the vacuum tube and silicon photomultipliers, the registration of the whole signal is not possible, and therefore sampling in the voltage domain using a predefined number of voltage levels is needed. The output signal is then recovered using ideas from the Tikhonov regularization (Tikhonov 1963, 1977) and compressive sensing (Candes et al2006, Donoho 2006) methods. The microchannel plates photomultipliers are the most promising in view of the application in the J-PET instrument due to the possibility of direct registration of a timestamp of each single photon. In the experimental study we will derive time resolutions of various configurations of the J-PET detector using different types of photomultipliers.

2 Materials and methods

In this work, we assume that the $\gamma$ photon interacts in the scintillator strip at time $\Theta$ and in the position $x$. We consider resolution for these reconstructions.

The time of the photon registration at the photomultiplier, referred to as $t_{r},$ is considered as a random variable, equal to the sum of three contributing values:

[TABLE]

where $t_{e}$ is the photon emission time, $t_{p}$ is the propagation time of the photon along the scintillator strip and $t_{d}$ is the photomultiplier transit time. Assuming that the times $t_{e},t_{p},t_{d},$ given in Eq. (1), are independent random variables with probability density functions (pdfs) denoted with $f_{t_{e}},f_{t_{p}},f_{t_{d}},$ respectively, the distribution function of $t_{r}$ is given as the convolution:

[TABLE]

In case of the ternary plastic scintillators used in the J-PET detector (Saint Gobain Crystals, Eljen Technology), the distribution of $t_{e}$ is well approximated by the following formula (Moszynski, Bergston 1977, 1979):

[TABLE]

where $\tau_{d}=1.5$ ns, $\tau_{r}=0.005$ ns and $\sigma_{e}=0.2$ ns, and $\kappa_{e}$ stands for the normalization constant. The values of the parameters $\tau_{d},\tau_{r},\sigma_{e}$ were adjusted in order to describe the properties of the light pulses from the BC-420 scintillator (Moskal et al2016, Saint Gobain Crystals). By definition in Eq. (2):

[TABLE]

Initial direction of flight of the photon in the scintillator is uniformly distributed. The photon on its way along the scintillator strip from the emission point to the photomultiplier may undergo many internal reflections whose number depends on the scintillator’s geometry and the photon’s emission angle. However, the space reflection symmetries of the cuboidal shapes, considered in this article, enables a significant simplification of the photon transport algorithm, without following photon propagation in a typical manner. The statistical modelling of this phenomena was presented in details in Ref. (Moskal et al2016) and the analytical function describing the distribution function $f_{t_{p}}$ may be expressed by the following formula:

[TABLE]

where $v$ is the speed of light in the scintillator strip, $\mu_{\text{eff}}$ is the effective absorption coefficient for the scintillator material and $\kappa_{p}$ the normalization constant. The $0\leqslant x\leqslant D$ is the longitudinal position of the emission point (see Fig. 2). The pdf function $f_{t_{p}}(t)$ in Eq. (4) is nonzero only for:

[TABLE]

where $t_{p}=\frac{x}{v}$ corresponds to the photon flying along the strip.

Finally, the time of registration $t_{r}$ is smeared using Gaussian distribution centered on the mean transition time $T_{d}$ and variance $\sigma_{d}^{2}$ estimated empirically:

[TABLE]

In this work, we assume that the signal registered at the photomultiplier output has the same functional dependence on the time as the $f_{t_{r}}$ function. We assume that the signal $y\in R^{N}$ is discretized by the oscilloscope. It is sampled in the constant time intervals denoted with $T_{s}.$ From the conditions Eq. (3) and (5), it follows that the registration time $t_{r}$ fulfils the inequality:

[TABLE]

It was assumed that the transition time $t_{d}\geqslant 0.$ Therefore, the $n^{th}$ time sample is given by:

[TABLE]

and the $n^{th}$ sample of the signal $y$ is given as:

[TABLE]

where $\beta(E,x)$ is a coefficient providing the scaling of the pdf function $f_{t_{r}}$ in order to obtain the voltage signal:

[TABLE]

The value of $\beta(E,x)$ depends on the energy deposited in the plastic scintillator during the $\gamma$ photon interaction ($\beta_{E}$ factor) and on the position of the $\gamma$ photon interaction along the strip ($\beta_{x}$ factor). The higher the value of deposited energy, the higher the value of $\beta_{E}$ parameter and higher the signal amplitude. The $\beta_{x}$ is necessary to describe absorption of photons propagating through the scintillator strip, since $f_{t_{p}}$ provides only information about the shape of the signal (see Eq. (4)). Hence, the closer to the left end of the scintillator, the smaller $x$ (see Fig. 2) and larger $\beta_{x}$. Contributions of $\beta_{E}$ to $\beta$ are the same for both ends of the strip but $\beta_{x}$ are different. Hereon, in order to simplify the notation of the parameter $\beta(E,x),$ we use only the symbol $\beta.$

2.1 Reconstruction of the interaction time and position

We denote the true values of time and position of $\gamma$ photon interaction with $\Theta^{0}$ and $x^{0},$ respectively, and the corresponding reconstructed values are denoted as $\hat{\Theta},\hat{x}.$ We add a random noise term $v_{(L,R)}$ to the signal $y_{(L,R)}$ at the left (L) and right (R) end of the strip. Hence a registered signals $\hat{y}_{(L)}$ and $\hat{y}_{(R)}$ may be expressed as:

[TABLE]

We assume that the noise $v_{(L)}$ and $v_{(R)}$ are uncorrelated and obey the same multivariate normal distribution:

[TABLE]

where $S$ is the covariance matrix of $\hat{y}_{(L)}$ and $\hat{y}_{(R)},$ and we introduce notation:

[TABLE]

According to the definitions of the theoretical ($y$) and registered ($\hat{y}$) signals, the reconstruction of $\hat{\Theta},\hat{x}$ may be pursued by minimization of the function:

[TABLE]

The solutions $\hat{\Theta},\hat{x}$ are found as:

[TABLE]

where hat denotes the estimators.

From Eqs. (9)-(10) and (12) it is seen that the error function $W$ is a positive-valued random variable. In order to determine $\Delta\hat{\Theta}$ we assume that error of time of interaction has normal distribution:

[TABLE]

where the $\sigma_{\Theta}$ is a searched time resolution of the J-PET instrument.

2.2 Determination of time resolution

In order to calculate the time resolution, $W$ has to be examined near the minimum, $(0,0).$ According to Eq. (12), the random variable $W(0,0)$ may be expressed as:

[TABLE]

The variance of $W$ in the minimum will be denoted hereafter as $\text{Var}[W_{\text{min}}].$ Using Eq. (11) and assuming the diagonality of matrix $S,$ yields:

[TABLE]

On the other hand, we may analyse the shape of the function $W$ in the two-dimensional space of time $(\Delta\hat{\Theta})$ and position $(\Delta\hat{x})$ errors. For the purpose of this work, we will consider only the $(\Delta\hat{\Theta})$ error, and therefore analyse $W$ in one dimension $(\Delta\hat{x}=0).$ Taylor series expansion of $W$ around $(0,0)$ is given as:

[TABLE]

It is evident that the first two coefficients ($\alpha_{0},\alpha_{1}$) are equal to zero and the quadratic approximation simplifies to:

[TABLE]

Under the assumption of normality of $\Delta\hat{\Theta}$ distribution, see Eq. (14), the random variable $W(\Delta\hat{\Theta},0)$ given in Eq. (18), has a $\chi^{2}$ distribution with the variance:

[TABLE]

The comparison of two formulas describing the $\text{Var}[W_{\text{min}}],$ in Eq. (19) and (16), enable us to determine time resolution, defined as the standard deviation $\sigma_{\Theta}:$

[TABLE]

2.3 Determination of coincidence resolving time

In order to facilitate the direct comparison with results published in the field of TOF-PET we will evaluate CRT based on the time resolution ($\sigma_{\Theta}$). In the first approximation CRT equals to $2.35\sqrt{2}\cdot\sigma_{\Theta}.$ However, a fundamental lower limit of the CRT is defined by the time spread due to the unknown depth-of-interaction (DOI) in a single scintillator. It should be stressed that this factor gains importance for large scintillator detectors as in e.g. J-PET. Since the interactions may occur with nearly equal probability along the whole thickness ($d$) of the plastic scintillator, time spread in a single scintillator may be well approximated by the uniform distribution with the width of $d/c,$ where $c$ denotes the $\gamma$ photon speed. This implies that the distribution of the time difference between two detected $\gamma$ photons has a triangle form with FWHM equal to $d/c.$ Therefore, the final value of CRT may be estimated with the formula:

[TABLE]

As seen from Eqs. (20)-(21), in order to evaluate $\sigma_{\Theta}$ and therefore CRT, one has to know the shape of pdf function $f_{t_{r}},$ to calculate the $\alpha_{2}$ coefficient, and also the errors of the signal registered on the photomultipliers, to calculate the covariance matrix $S$. Determination of the shape of $f_{t_{r}}$ was discussed in the previous section. In the next section we will analyse the sources of errors in the signals $\hat{y}_{(L)},\hat{y}_{(R)}.$

2.4 Analysis of registered signals errors

The noise contribution to the signals registered on the left ($\hat{y}_{(L)}$) and right ($\hat{y}_{(R)}$) side of the scintillator strip is the same, and therefore in this section we will skip the $L,R$ indices. In further analysis we assume that the noise signal $v,$ see Eq. (9), is defined as a sum of two components:

[TABLE]

where $v_{p}$ describes the perturbations of the pdf function $f_{t_{r}},$ based on limited number of input photon signals, and $v_{r}$ stands for the signal recovery noise. The latter component is introduced by the procedure of signal recovery based on the limited number of registered samples of the signal in the voltage domain. The problem of signal recovery was widely discussed in Ref. (Raczynski et al2015a, b). We assume that the noises $v_{p}$ and $v_{r}$ are uncorrelated and normally distributed with covariance matrices $S_{p}$ and $S_{r},$ respectively. Thus, one may write that:

[TABLE]

The exact values of $v_{p}$ and $v_{r}$ depend on the type of the photomultiplier applied. In this work we consider two types of photomultipliers:

• •

PMT - vacuum tube photomultiplier treated as a basic one in the current J-PET prototype,

• •

MCP - microchannel plates photomultiplier.

It should be underlined that the following analysis does not include the silicon photomultipliers. We have provided the extensive research of the possibilities of the application of the silicon photomultiplier in the J-PET tomograph in our previous study in Ref. (Moskal et al2016).

Noises $v_{p}$ and $v_{r}$ are mainly influenced by the width of the single photoelectron contributing to the final output signal, and quantum efficiency of the photomultiplier. The most distinctive feature of the MCP photomultiplier is the capability of the registration of arrival time of each photon. Thus, the output signal may be evaluated by using a model of the single photon. For all types of photomultipliers we use the Gaussian model (Bednarski et al2014) for shape of signal of single photoelectron, with the width $\sigma_{p}.$ In the experimental section we will optimize $\sigma_{p}$ parameter for the MCP photomultiplier, aiming to minimize the $v_{p}$ noise. The quantum efficiency may be used directly to estimate the number of photoelectrons induced in the photomultiplier, $N_{p}.$ In the following we will apply $N_{p},$ to model the total output signal.

It is worth noting that $v_{r}$ vanishes in the case of MCP photomultiplier. There is no need to recover the output signal since all arrival times of photons are registered. In the following we will shortly describe the noises $v_{p}$ and $v_{r}.$

2.4.1 Analysis of $v_{p}$.

The registered signal $y$ affected only by the $v_{p}$ noise will be denoted with:

[TABLE]

The signal $\tilde{y}$ consists of $N_{p}$ signals from individual photoelectrons:

[TABLE]

As mentioned in Section 2.4, signal from single photoelectron $\tilde{y_{k}}$ is assumed to be a Gaussian function:

[TABLE]

where $t_{r}^{k}$ is a random variable with $f_{t_{r}}$ distribution, that denotes the $k^{th}$ photon’s registration time.

We aim to calculate the diagonal elements of the covariance matrix $S_{p}:$

[TABLE]

where

[TABLE]

According to the Eq. (24):

[TABLE]

Estimates of the $\text{Var}((\tilde{y}(n))$ and $\text{Bias}(\tilde{y}(n))$ were introduced in Refs. (Rosenblat 1956, Simonoff 1996). Assuming that the underlying pdf function $f_{t_{r}}$ is sufficiently smooth, and that $\sigma_{p}\rightarrow 0$ with $N_{p}\sigma_{p}\rightarrow\infty$ as $N_{p}\rightarrow\infty,$ the Taylor series expansion gives:

[TABLE]

where $f_{tr}^{\prime\prime}(t^{(n)})$ is a second derivative of the pdf function $f_{tr}(t^{(n)}).$ Above approximations may be inaccurate for finite $N_{p}.$ The number of registered photoelectrons $N_{p}$ is of the order of hundreds, and the detailed discussion is given in Sec. 3.1. Therefore, a new method to evaluate the $\text{Var}((\tilde{y}(n))$ and $\text{Bias}(\tilde{y}(n))$ for finite $N_{p}$ should be proposed. During this study the novel concept of the estimation of a requested statistics has been developed. The method has been described in great details in the Appendix and it was shown that the values of $\text{Var}(\tilde{y}),\text{Bias}(\tilde{y})$ may be estimated as:

[TABLE]

where $\lambda$ is the parameter defining the range of the second argument of function $\Phi$:

[TABLE]

and $F_{t_{r}}(t^{(n)})$ is the cumulative distribution function of $f_{t_{r}}(t^{(n)})$ calculated at $t^{(n)}.$ Discussion of formulas (32, 33) is given in the Appendix.

It should be underlined that both estimation methods, proposed (Eqs. (32, 33)) and based on Taylor series approximation (Eqs. (30, 31)), have the same asymptotic properties. It may be shown that for $\sigma_{p}\rightarrow 0$ with $N_{p}\sigma_{p}\rightarrow\infty$ as $N_{p}\rightarrow\infty:$

[TABLE]

2.4.2 Analysis of $v_{r}$.

Denote the signal $y$ affected only by the $v_{r}$ noise as:

[TABLE]

The recovery process takes place only provided the complete output $\tilde{y}$ is registered on a photomutliplier. If times of photon’s arrival are registered, as in the MCP photomultiplier, the $v_{r}=0.$ Recovery of the signal $\hat{y}$ is carried out only for the PMT photomultiplier.

The details of signal recovery process were given in Ref. (Raczynski et al2015a), and here only the main points will be recalled. The evaluation of the signal $\hat{y}$ requires two steps: (i) recovery of the sparse expansion $\hat{x}$ and (ii) calculation of $\hat{y}$ based on the $\hat{x}$. The relation between the solution $\hat{y}$ and its sparse representation $\hat{x}$ is linear:

[TABLE]

where $A$ is an orthonormal matrix. As it was shown in Ref. (Raczynski et al2015a), from the Bayes theory the properties of regularized solution can be found, in particular its covariance matrix, denoted hereafter as $S_{r(x)}$, may be easily derived:

[TABLE]

where $P$ is the covariance matrix of the sparse signals $x$, and $M$ denotes the number of registered samples of the signals $y,\sigma$ is the standard deviation of the measurement error. Finally, based on Eq. (36), the covariance matrix $S_{r}$ is given:

[TABLE]

3 Experimental results

3.1 Experimental setup

In this section we investigate the accuracy of the proposed method for evaluation of the time resolution and CRT. The model is validated by performing the experiment with a single detection module of the J-PET scanner built out from the BC-420 plastic scintillator strip, with dimensions of 5 x 19 x 300 mm, read out at two ends by the Hamamatsu R4998 (PMT) photomultipliers. Our experimental setup is depicted in Fig. 3.

Measurements are performed using $\gamma$ photon from the 22Na source placed inside the lead collimator between the scintillator strip and the reference detector. The reference detector consists of a small scintillator strip with a thickness of 4 mm. Collimated beam emerging through 1.5 mm wide and 20 cm long slit is used for irradiating desired points across the strip. In order to detect the event, a coincident registration of signals from the PM1 and a reference detector is required. Such trigger conditions enable us to select precisely the annihilation quanta reducing the background from the deexcitation photon (1.27 MeV) to the negligible level (Moskal et al2014b). The time of triggering by the reference detector is used to estimate the event arrival time. The constant electronic time delay between the true event time and measured arrival time to the reference detector does not influence the time resolution and is shifted to zero. The full waveforms of PMT signals are sampled using the Lecroy SDA 6000A oscilloscope running at a 20 GSps sampling rate.

In our previous studies it was shown that the time resolution is fairly independent of the irradiation position (Moskal et al2014b). Therefore, we determine the time resolution and CRT of the J-PET scanner in one position, at the center of the strip ($x=$ 15 cm). In order to evaluate the experimental value of time resolution and CRT a data set of $10^{4}$ pairs of signals from PM1 and PM2 registered in coincidence was analyzed. In the first step, for each pair of fully sampled signals from the left and right ends of the strip, $\tilde{y}_{(L)}$ and $\tilde{y}_{(R)},$ a front-end electronic device probing signals at four voltage levels, both at the rising and falling slope, was simulated. The signals $\hat{y}_{(L)}$ and $\hat{y}_{(R)}$ were recovered using 8 samples of signals $\tilde{y}_{(L)}$ and $\tilde{y}_{(R)}$ registered by an oscilloscope, according to the method descried in Sec. 2.4.2. For the $k$-th pair of the recovered signals $\hat{y}_{(L)}$ and $\hat{y}_{(R)},$ the energy of an event may be estimated based on arithmetic mean of a number of photoelectrons registered at the left and right sides of the scintillator (Moskal et al2014b) and is proportional to the sum of integrals of recovered signals $\hat{y}_{(L)}$ and $\hat{y}_{(R)}.$ On the other hand, for the $k$-th pair of $\hat{y}_{(L)}$ and $\hat{y}_{(R)},$ the reconstruction of time ($\hat{\Theta}_{k}$) and position ($\hat{x}_{k}$) was pursued by minimization of the function $W$ in Eq. (12). The value of $\sigma_{\Theta}$ was calculated as the standard deviation of the empirical distribution of $\hat{\Theta}_{k}$ and was equal to about 80 ps. The corresponding value of CRT calculated based on Eq. (21), for a scintillator strip with a thickness of 19 mm, was equal to 275 ps. This value of CRT will be treated as the reference for the proposed approach. For clarity of the presentation, we will calculate in Sec. 3 only the CRT parameter.

According to the scheme presented in Sec. 2.4, the evaluation of the time resolution and CRT of the PET system requires the investigation of the parameter $\alpha_{2}$ and covariance matrix $S$ (see Eq. (20)). The values of these parameters vary for a different type of applied photomultipliers and are also sensitive to the position of the point of $\gamma$ photon interaction along the scintillator strip. The values of the parameters will be provided in Sec. 3.1 and 3.2.

In order to model the signal at the photomultiplier’s output, the parameters of three pdf functions $f_{t_{e}},f_{t_{p}}$ and $f_{t_{d}},$ defined in Eq. (2), (4) and (6), respectively, must be known. It is worth noting that only the last pdf function, $f_{t_{d}},$ describes the unique properties of a given type of the photomultiplier. The $\sigma_{d}$ is delivered by the photomultiplier’s producer, for Hamamatsu R4998 photomultiplier (PMT) $\sigma_{d}=$ 68 ps and for MCP photomultiplier $\sigma_{d}=$ 40 ps (Hamamatsu 2016). However, as our initial tests show, there is a negligible influence of the $\sigma_{d}$ value on the performance of the reconstruction method of $\gamma$ photon interaction moment.

The measurement provides a discrete signal $y$ at the photomultiplier. Repeating the measurement of its time under the same condition yields a set of the acquisition times of photoelectrons, cf. Eq. (24). The signal $y$ consists of $N_{p}$ Gaussian shaped signals of single photoelectrons. In Fig. 4 an example of the single photoelectron signal registered with PMT photomultiplier and its Gaussian fit are shown. The signals are marked with blue and red curves, respectively. The standard deviation $\sigma_{p}$ of this function is reported in Ref. (Bednarski et al2014) to be equal to 300 ps in the case of PMT photomultiplier.

However, a different number of photoelectrons ($N_{p}$) are registered due to the different quantum efficiencies. In the following we will shortly recall the main results of our earlier works enabling us to estimate properly the number $N_{p}.$ Light yield of plastic scintillators amounts to about 10000 photons per 1 MeV of deposited energy. The 511 keV $\gamma$ photon may deposit maximally 341 keV via Compton scattering (Szymanski et al2014), which corresponds to the emission of about 3410 photons. On the other hand in order to decrease the noise due to the scattering of $\gamma$ photon inside patient’s body, the minimum energy deposition of about 200 keV is required (Moskal et al2012). Therefore, the range of the number of emitted photons discussed hereafter in this article amounts 2000 to 3410. The experiments conducted with PMT photomultipliers revealed that about 280 photoelectrons are produced from the emission of 3410 photons (Moskal et al2014). According to the preselected range, from 2000 to 3410 photons, the average number of emitted photons is about 2700. This number corresponds to $N_{p}=$ 220 registered photoelectrons equipped with the PMT photomultiplier. Since the CRT of the J-PET system will be determined at the center of the strip, the numbers of photoelectrons $N_{p}$ contributing to the signals induced on the left and right scintillator ends are the same and are equal to 110.

As mentioned at the beginning of Sec. 2, the values of $\tau_{d},\tau_{r},\sigma_{e}$ of the $f_{t_{e}}$ pdf function were adjusted based on the experimental studies with a single BC-420 scintillator strip. We have provided numerous tests for various strips of the BC-420 scintillator type and we found that the values of the estimated parameters of the $f_{t_{e}}$ pdf function were consistent within the measurement errors. Therefore, the signals evaluated based on the proposed model, presented in Sec. 2.2, have shapes very similar to those registered during the experiment via oscilloscope (see Fig. 5). In Fig. 5, the theoretical signal $y$ at the center of the strip, evaluated from Eq. (8), is presented. The parameter $\beta$ (see formula (8)) was selected in such way that the amplitude of the signal is equal to the mean amplitude of signals registered at the center of the strip ($x=$ 15 cm). The analytical solution for $f_{t_{r}}$ function is difficult to find due to the internal convolution in $f_{t_{e}}$ function (see Eq. (2)). Therefore, the numerical evaluation of a convolution operation was applied. The signals $y$ and $\tilde{y}$ in Fig. 5 are shown in the discrete domain for discrete time samples and the curves connecting points are plotted to guide ones eye.

The information about the signal $y$ may be directly applied to evaluate the value of the parameter $\alpha_{2}.$ In this work we are interested only in determination of CRT of the J-PET system and we assume that the position of the $\gamma$ photon interaction is known exactly (see Eq. (18)). Therefore, for a fixed position of the interaction, the signal $y$ may be shifted only in time domain due to the error of time measurement $\Delta\Theta.$ For $\Delta\Theta=$ 0 the theoretical and registered signals overlap and $W(\Delta\Theta,0)=$ 0, see Eq. (12). In order to evaluate $\alpha_{2},$ the error $\Delta\Theta$ was varied from -1 to 1 ns. For each value of $\Delta\Theta,$ function $W(\Delta\Theta,0)$ was evaluated based on the shape of signal $y$ shown in Fig. 5. The resulting, experimental function $W(\Delta\Theta,0)$ is presented in Fig. 6 with blue curve (see also Eq. (12)).

According to Eq. (18), the experimental function $W(\Delta\Theta,0)$ may be approximated near $\Delta\Theta=$ 0 with the quadratic function. The quadratic approximation of the $W(\Delta\Theta,0)$ function is marked in Fig. 6 with the red curve and the coefficient of the second order polynomial function is equal to 11.2 $\frac{\text{V}^{2}}{\text{ns}^{2}}.$

3.2 Verification of signal $\tilde{y}$ estimation method

According to the assumptions in Sec. 2.4, two main contributors to the signal’s noise are $v_{p}$ and $v_{r}.$ The $v_{r}$ was estimated in Ref. (Raczynski et al2015a) and will be recalled at the end of this point. In this section, a detailed study of the approximation method of the $\text{Var}(\tilde{y}(n))$ and $\text{Bias}(\tilde{y}(n))$ will be carried out. The proposed method, see Eq. (32) and (33), will be compared with the well known approximation technique based on the Taylor series expansion, see Eq. (30) and (31). As the reference for the results of both analytical approaches, the Monte-Carlo (MC) simulation will be provided.

The MC simulation was carried out for the constant number of photoelectrons $N_{p}=$ 220, registered by the PMT photomultiplier. In order to simulate the $\text{Var}(\tilde{y}(n))$ and $\text{Bias}(\tilde{y}(n)),$ only one timestamp of the original signal $y,$ corresponding to the maximum value of 0.6 V (see Fig. 5), was used. The analysis of the maximum value in signal $y$ allows one to evaluate the main contribution in the covariance matrix $S_{p};$ as seen from Fig. 9, the location of maximum value of the signal $y$ corresponds to the location of maximum value on diagonal of the covariance matrix $S_{p}$. The maximum value of the original signal $y$ is observed in the sample $n=$ 60 (see Fig. 5). In the first step of MC simulation the random values of photons registration times $t_{r}^{k}$ ($k=1,2,...,N_{p}$) were selected according to the $f_{t_{r}}$ distribution. Next, the values of all $N_{p}$ functions $\tilde{y_{k}}(60)$ were evaluated based on the Eq. (25) and summed up giving $\tilde{y}(60).$ The above-mentioned procedure was repeated $10^{6}$ times for different values of $\sigma_{p}$ from 50 ps to 750 ps with step 25 ps. The range of $\sigma_{p}$ has been selected after a preliminary calculations taking into account the expected number of registered photoelectrons in the J-PET scenario. Based on the large number of samples of $\tilde{y}(60),$ the accurate estimation of bias and variance was possible. The resulting $\text{Bias}^{2}(\tilde{y}(60))$ and $\text{Var}((\tilde{y}(60))$ are shown in Fig. 7 and 8, respectively.

The reference values of $\text{Bias}^{2}(\tilde{y}(60))$ and $\text{Var}(\tilde{y}(60)),$ obtained with MC simulation, are marked with red curves in Fig. 7 and 8, respectively. An approximation of $\text{Var}(\tilde{y})$ for the proposed method and method based on the Taylor series expansion (blue and green curves, respectively) are very similar to reference curve for small values and tend to differ for larger values of $\sigma_{p}.$ However, in the most interesting region, for $\sigma_{p}$ equal to about 300 ps, the proposed method is more accurate than the Taylor series based method and the values of $\text{Var}(\tilde{y})$ are equal to 6.5$\times 10^{-3}$ V2 and 7.0$\times 10^{-3}$ V${}^{2},$ respectively (the reference value of $\text{Var}(\tilde{y})$ from MC simulation is equal to 5.1$\times 10^{-3}$ V2). Comparison of the $\text{Bias}^{2}(\tilde{y})$ and $\text{Var}(\tilde{y})$ curves reveals the fundamental relation between variance and bias. The variance dominates for smaller values of $\sigma_{p}$ and becomes comparable with bias for $\sigma_{p}$ at the level of about 500 ps (compare two reference, red curves in Fig. 7 and 8). For $\sigma_{p}$ larger than 500 ps, the total error is mostly influenced by the bias. It is worth noting that in that case the Taylor series based method significantly underestimates the values of $\text{Bias}^{2}(\tilde{y}),$ see Fig. 7, which leads to the underestimation of the overall error.

3.3 Evaluation of time resolution of the J-PET system

In the first step we compare the covariance matrices $S_{p}$ and $S_{r}$ (see Eqs. (22) and (23)) according to the description in Sec. 3.1 and using our previous study (Raczynski et al2015a). The resulting values of the diagonal elements of $S_{p}$ and $S_{r}$ are shown in Fig. 9.

Theoretical values of $S_{p}$ were evaluated as in Sec. 3.2 for PMT and the results are marked in Fig. 9 with blue curve. The values of $S_{r}$ elements are marked with green curve in Fig. 9. The comparison of the resulting characteristics with the shape of the pdf function $f_{t_{r}},$ presented in Fig. 5, indicates that the reconstructed errors $S_{p}$ and $S_{r}$ are highly related to the signal value. The maximal values of the diagonal elements of $S_{p}$ and $S_{r}$ occur near to the maximum of signal $y$ (Fig. 5). On the other hand, the analysis of the characteristics plotted in Fig. 9 reveals that the error introduced by the limited number of photoelectrons in the registered signal ($S_{p}$) is a dominating factor.

In order to compare the reconstructed values of the covariance matrices $S_{p}$ and $S_{r},$ we use the trace (Tr) of the covariance matrix since the diagonality is assumed. The values of $\textnormal{Tr}(S_{p})$ and $\textnormal{Tr}(S_{r})$ for different photomultipliers type are gathered in Tab. 1. In the following we will analyse the value of $\textnormal{Tr}(S_{p})$ as the function of the number of registered photoelectrons ($N_{p}$) and standard deviation of the single photoelectron signal ($\sigma_{p}$). The resulting characteristics of $\textnormal{Tr}(S_{p})$ as a function of $\sigma_{p}$ are shown in Fig. 10.

The values of $\textnormal{Tr}(S_{p})$ were calculated for three specified numbers of registered photoelectrons 220, 350 and 700 and are marked in Fig. 10 with blue, red and black curves, respectively. The smallest value of $N_{p}$ is specific for PMT photomultiplier, as mentioned in Sec. 3.1. The highest number $N_{p}=$ 700 indicates the maximal number of registered photoelectrons in the experimental scenario, and was selected in order to demonstrate the best theoretical resolution of the J-PET. Results in Fig. 10 show that all the $\textnormal{Tr}(S_{p})$ functions, evaluated for a given number of registered photoelectrons, have a minimum. The shape of the $\textnormal{Tr}(S_{p})$ functions illustrates the fundamental trade-off between variance and bias, as mentioned in Sec. 3.2. Hence, for given number of registered photoelectrons, it is possible to adjust the optimal value of $\sigma_{p}$, denoted hereafter with $\sigma_{p(opt)}.$ Comparison of the $\sigma_{p(opt)}$ values, for three $N_{p}$ numbers in Fig. 10, shows that the larger the number of registered photoelectrons the smaller the value of $\sigma_{p(opt)}.$ For instance, for the PMT photomultiplier that registers 220 photoelectrons on average, the minimum error occurs for the $\sigma_{p(opt)}=$ 500 ps (blue curve in Fig. 10). In the case of PMT photomultiplier the $\sigma_{p}$ is not a variable and has fixed value of about 300 ps. However, the MCP photomultiplier registers timestamps of the signal instead of the complete signal. Therefore, the value of $\sigma_{p(opt)}$ of each contributing signal may be adjusted accordingly to the number of registered timestamps ($N_{p}$). In that sense, the optimization of the $\sigma_{p(opt)}$ value for MCP photomultiplier may be provided. Simulations using $f_{t_{r}}$ function provide $N_{p}^{-4.1}$ dependence of $\sigma_{p(opt)}$ on the number of photoelectrons.

In general, the MCP photomultiplier is capable to register all the timestamps of the photons reaching the scintillator end. In order to account for possible inefficiency of the MCP, we determine the characteristics of the J-PET equipped with the MCP in the range $100\leqslant N_{p}\leqslant 700.$ First, for a given number $N_{p},$ the optimal value of $\sigma_{p(opt)}$ was estimated based on the characteristics of $\textnormal{Tr}(S_{p})$ (see Fig. 10). Next, the matrix $S_{p}$ was calculated based on the proposed technique, see Eq. (32) and (33). Finally, $\sigma_{\Theta}$ was evaluated based on Eq. (20). In the case of MCP, $\textnormal{Tr}(S_{r})=$ 0, since the output signal is given directly based on the measured timestamps and assumed shape of the single photoelectron signal. Resulting characteristics of CRT is given with red solid line in Fig. 11. The presented values of CRT take into account an additional smearing of the time due to the unknown depth of interaction in a scintillator strip with a thickness of 19 mm, see Eq. (21) for details.

Additionally, the CRT calculated for $\sigma_{p}=$ 300 ps (a shape of single photoelectron signal characteristic for PMT photomultiplier), including also the $v_{r}$ error introduced by the signal recovery procedure, $\textnormal{Tr}(S_{r})=$ 0.22 V2 reported in Ref. (Raczynski et al2015a), is marked with blue dashed line in Fig. 11. The theoretical value of CRT for the PMT photomultiplier is marked with full circle on the blue curve, for $N_{p}=$ 220. The theoretical CRT of the J-PET scanner with the PMT photomultiplier is about 290 ps and agrees with the experimental value of CRT, reported to be about 275 ps (Raczynski et al2015a). For fixed value of the quantum efficiency (equivalent to the number $N_{p}$), further improvement of CRT is possible by the application of the MCP photomultipliers. In the wide range of numbers of registered photoelectrons shown in Fig. 11, improvement of about 30 ps is observed (red and blue dashed curves). The presented results show that the best theoretical CRT of the J-PET scanner, with 30 cm long strips, estimated for the MCP photomultiplier that is capable to register all timestamps of arrival for 700 photons, is at the level of 170 ps. The main results of the study as well as the parameters of the analysed photomultipliers, are summarized in Tab. 1.

4 Extension of the proposed method for conventional PET systems

The proposed framework for calculation of the time resolution and CRT may also be applied for the state of the art PET scanners equipped with crystal scintillators. To this purpose, in order to estimate the time resolution, function $W,$ introduced in Eq. (12), has to be adopted to the new situation. First of all, during the reconstruction of the time of $\gamma$ photon interaction ($\hat{\Theta}$), only signals from one side of the crystals are acquired; in J-PET scanner two signals are registered at both ends of scintillator. Moreover, due to the small size of the crystals, only reconstruction of the interaction time $\hat{\Theta}$ is being carried out; function $W$ is one-dimensional. Therefore, the function $W$ is defined as:

[TABLE]

and finally, the standard deviation $\sigma_{\Theta}$ is given with the formula:

[TABLE]

which differs from the Eq. (20) only with a factor $\sqrt[4]{2}.$ In order to evaluate time resolution and CRT, one has to evaluate the distribution time of the photon registration at the photomultiplier ($t_{r}$), to calculate the $\alpha_{2}$ coefficient, and also the errors of the signal registered on the photomultipliers, to calculate the covariance matrix $S.$

In case of PET scanners with inorganic crystal scintillators, the description of the photon registration time $t_{r}$ at the photomultiplier includes two random components $t_{e}$ and $t_{d}:$

[TABLE]

In comparison to the J-PET, here the propagation time of the photon along the scintillator ($t_{p}$) may be neglected due to small size of the single crystal. Therefore, the only difference in calculation of $t_{r}$ is the evaluation of distribution of $t_{e}.$ The main parameter that governs the speed of light emission after the absorption of a $\gamma$ photon is the decay time. Crystal materials show decay pulse shapes that are single- or multi-exponential. For example, for the BGO crystals the bi-exponential shape of the distribution of time $t_{e}$ is observed (Seifert et al2012).

On the other hand, determination of time resolution, defined in Eq. (39), requires the information about the covariance matrix $S$ of the registered signal $\hat{y}.$ In Sec. 2.4 we derived the analytical description of the main components of the covariance matrix $S,$ i.e. matrices $S_{p}$ and $S_{r}.$

The formula for calculation of elements of $S_{p}$ matrix, describing the perturbations of the distribution function $f_{t_{r}}$ based on a limited number of registered photoelectrons, given in Eq. (26), may be also applied to the PET scanners with crystal scintillators. The only differences are in the parameters describing the shape of the $f_{t_{r}}$ distribution function and the expected number of photoelectrons ($N_{p}$) while including the light yield of crystal scintillators.

The latter component, matrix $S_{r},$ is introduced by the procedure of signal recovery based on the limited number of registered samples of the signal in the voltage domain. The J-PET system involves a four-threshold sampling method to generate samples of a signal waveform. An example of a similar electronic system for probing the signals in a voltage domain, coupled with experimental setup equipped with LSO crystals was developed in Ref. (Kim et al2009). The waveforms of signals were read-out by the oscilloscope and the electronic system for probing these signals in a voltage domain with four thresholds was applied to reconstruct the pulse shape. This scheme allows to evaluate all the parameters required to calculate the signal recovery error, according to the formula given in Eq. (38), for the PET system with crystal scintillators.

5 Conclusions

In this paper, we introduced a new method of estimation of time resolution and CRT of the J-PET system using only simulations which were tested based on the data from a single detector module. This is particularly useful for design of expensive device. In case of J-PET tomograph the most expensive part of the system are the photomultipliers. In this work two types of photomultipliers were simulated: the vacuum tube photomultipliers and microchannel plates.

The basic idea of the method is the use of the statistical nature of the whole signal acquisition process. We have highlighted three statistical phenomena: the emission of photons in the scintillator strip, the propagation of light pulses along the strip and registration of light in photomultipliers. Parameters of the probability density functions were selected in order to describe properly light pulses from the plastic scintillator BC-420.

An important aspect of our work concerns the statistical analysis of an error of reconstruction of the probability density function based on the set of single photoelectrons signals. In this work dependences of an overall variance and bias on the number and width of the single photoelectron signals were evaluated. The proposed estimation method was validated by using the Monte Carlo simulation and it was shown that obtained results are consistent. Moreover, the proposed technique was demonstrated to be more accurate than the approach from the literature (Rosenblat 1956, Simonoff 1996). The developed estimation scheme is general and may be incorporated elsewhere.

In the experimental section, the method of time resolution and CRT estimation was tested using signals registered by means of the single detection module of the J-PET scanner. In order to evaluate a CRT of the J-PET detector, we have incorporated the method described in Ref. (Raczynski et al2015a). In the cited work, the CRT obtained with the experimental scheme with vacuum tube photomultipliers, was reported to be equal to about 275 ps. Our calculation shows that the application of the proposed estimation method can give very similar result of about 290 ps. The consistency of the experimental and theoretical results, obtained for the J-PET scanner equipped with the vacuum tube photomultipliers suggests that the estimated CRTs for other photomultipliers are reliable. The determined CRTs for the detector with microchannel plates amount to 215 ps and 170 ps assuming the 50$\%$ and 100$\%$ quantum efficiency of photomultiplier, respectively.

Future work will address investigation of other aspects of signal acquisition process by using the proposed statistical model, e.g. the influence of the parameters of distribution of the photon emission time on time resolution. In this study, the parameters of the distribution were selected in order to describe the properties of light signals observed in the BC-420 plastic scintillator. However, our group develops a novel type of plastic scintillator and examines the influence of the chemical composition of the plastic scintillator on the overall performance of the J-PET detector (Wieczorek et al2015a, 2015b, 2016). Application of the proposed model to that task, enables us to use information about the shape of the distribution of the time of photon emission directly to predict the CRT of the J-PET detector.

We acknowledge technical and administrative support of A. Heczko, M. Kajetanowicz, W. Migdał and the financial support by the Polish National Center for Development and Research through grant INNOTECH-K1/IN1/64/159174/NCBR/12, the Foundation for Polish Science through MPD programme, the EU and MSHE Grant No. POIG.02.03.00-161 00-013/09, Doctus - the Lesser Poland PhD Scholarship Fund, and Marian Smoluchowski Kraków Research Consortium ”Matter-Energy-Future”.

Appendix A Kernel density estimation

The function $\tilde{y_{k}},$ describing the $k^{th}$ signal from a single photoelectron, given in Eq. (25), may be approximated with:

[TABLE]

where $n=1,2,...,N,$ and $\lambda$ contributes to the signal width. The probability that the random variable $\tilde{y_{k}}(n)$ is equal to the specified value, may be calculated based on the previously introduced function $\Phi$, see Eq. (34). In particular, the probability that the random variable $\tilde{y_{k}}(n)=0$ is equal to $1-\Phi(t^{(n)},\lambda\sigma_{p});$ the $k^{th}$ registration time $t_{r}^{k}$ is out of range $(t^{(n)}-\lambda\sigma_{p},t^{(n)}+\lambda\sigma_{p}),$ see the second case in Eq. (40). Denoting the first case in Eq. (40) with $u_{k}:$

[TABLE]

we may write that for $n=1,2,...,N,$ the expected value of $\tilde{y_{k}}(n)$ is equal:

[TABLE]

and the variance of $\tilde{y_{k}}(n)$ is equal:

[TABLE]

In order to simplify the further calculations the following assumption is proposed. Note that in most interesting cases the range $(t^{(n)}-\lambda\sigma_{p},t^{(n)}+\lambda\sigma_{p}),$ is narrow in comparison to the estimated pdf function $f_{t_{r}}$ domain. Therefore, the pdf function $f_{t_{r}}$ is considered to be uniform in the range $(t^{(n)}-\lambda\sigma_{p},t^{(n)}+\lambda\sigma_{p}):$

[TABLE]

It is worth noting that the smaller is the ratio of the single- to overall signal width, the better is the performance of the proposed approximation method.

Under the assumption in Eq. (44), required moments in Eqs. (42) and (43), $E[u_{k}(n)],$ $E[u_{k}(n)]^{2}$ and $\text{Var}(\tilde{u_{k}}(n))$, can be easily derived. After some simple calculations the equations for the expected value and the variance of the random variable $\tilde{y_{k}}(n)$ are given by formulas:

[TABLE]