Low-Pass Filtering Method for Poisson Data Time Series

: Problems of digital processing of Poisson-distributed data time series from various counters of radiation particles, photons, slow neutrons etc. are relevant for experimental physics and measuring technology. A low-pass ﬁltering method for normalized Poisson-distributed data time series is proposed. A digital quasi-Gaussian ﬁlter is designed, with a ﬁnite impulse response and non-negative weights. The quasi-Gaussian ﬁlter synthesis is implemented using the technology of stochastic global minimization and modiﬁcation of the annealing simulation algorithm. The results of testing the ﬁltering method and the quasi-Gaussian ﬁlter on model and experimental normalized Poisson data from the URAGAN muon hodoscope, that have conﬁrmed their effectiveness, are presented.


Introduction
The article proposes a low-pass filtering method for Poisson-distributed data time series, based on a specially developed digital low-pass filter with finite impulse response (FIR filter), with gain equal to one at zero frequencies and non-negative weighting factors.
Here, low-pass filtering is applied in order to reduce noise in Poisson-distributed data to ensure the recognition of emerging fluctuations of mathematical expectations in them. Poisson-distributed, or Poisson data are found in various physical systems, for example, related to the heliosphere and magnetosphere of the Earth; the fluctuations of mathematical expectations of these data may contain information regarding the structures and characteristics of these systems.
A particular feature of the Poisson data origin is that they contain sufficient noises; it is known, for example, from [1] that their variance is numerically equal to mathematical expectation. Noise reduction in Poisson data can be achieved using common FIR filters [2,3], to which, within the framework of this article, we refer the filters based on commonly used windowing techniques, frequency sampling and inverse Fourier transforms [4,5]. However, there are a number of scientific and technical problems for which their application is not fully effective, for example, (1) recognition of small (in size and duration) mathematical expectation fluctuations in Poisson datasets; (2) digital processing of Poisson data with small mathematical expectation values.
Common FIR filters can potentially be used for the mentioned tasks, and their synthesis can be implemented according to given dimensions and cutoff frequencies. The synthesis procedures for common FIR filters are, in essence, the variants of approximation procedures for the specified species frequency response (FR) types; the accuracy of the FR approximations depends on the specified dimensions for the synthesized filters. Obviously, at large dimensions, the accuracy of these approximations is high and the errors in the resulting cutoff frequencies are small. For the case of small dimensions, the approximation accuracy turns out to be low and, as a consequence, cutoff frequencies are realized with significant errors which prevent low-pass filtering. We can assume that the filtering procedure proposed here should be performed by filters with low dimensions and cutoff frequencies and with gain values equal to one in order to avoid mathematical expectation distortions, and with non-negative weight factors in order to provide non-negativity of filtering results taking into account the Poisson property of the data.
The indicated problem leads to the need to formulate the synthesis problem for a special digital low-pass FIR filter, which takes into account the requirements-restrictions on dimensionality, cutoff frequency, gain at zero frequencies, and weighting factors.
Here, a FIR filter is proposed, which is further denoted as a quasi-Gaussian filter, the frequency response of which is formed on the basis of approximating a Gaussian function and ensuring the implementation of the mentioned constraints conditions using a special optimization method.
Gaussian filters, the frequency response of which is implemented based on the approximation of the Gaussian function, are widely used in modern scientific and technical practice [6,7]. However, as a rule, the known variants of Gaussian filters with the approximation of the frequency response do not take into account the above-mentioned conditions (restrictions).
Problems of digital processing of Poisson data time series from muon counters in muon detectors and telescopes [8], counters of elementary particles of alpha-beta-gamma radiation, photon counters, slow neutrons, etc. [9], taking into account their specificity, are relevant for experimental physics. Digital processing of Poisson data, including the Gaussian filtering application, can be outside of experimental physics, for example, in medical technology for imaging blood vessels and tumor therapy with particle beams, in measuring technology for tribological studies of the surfaces of metal parts, in astronomy for gamma telescopes, in muon tomography for recognizing cavities in rocks, and building structures and many other applications.
One of the applications of the designed filter proposed here is the digital processing of the data from the URAGAN muon hodoscope (MH) designed by NRNU MEPhI [10,11]. The MH is a computerized measuring device that estimates the intensities of muon fluxes by counting the number of elementary particles-muons-registered in its detector for a set of solid angles with a set time step. Within the framework of this article, MH can be interpreted as a distributed set of muon counters, consisting of primary and secondary information converters.
From each primary MH transducer, the initial Poisson data-time series of random non-negative integers N(Tk)-the quantities of Poisson-distributed events recorded in a given solid angle at time intervals (T(k − 1), T(k − 1) + T 0k ), k = 1, 2, . . . , k 0 , where T = 1 minute. Due to the features of the MH design, registration intervals T 0k are random with a uniform distribution law in the range T 0 min ≤ T 0k ≤ T 0 max < T.
From each secondary MH transducer, the 1-minute-sampled normalized Poisson data Y(Tk) are generated for a given solid angle by reducing to one second and calculating the averaged normalized Poisson data Y(T 0 n) with an hourly discreteness according to the following relations:  One-dimensional FIR filter synthesized here is built according to the following difference equation: where r 0 = s 0 + 1 is the FIR filter dimension, a T = (a 0 , a 1 , . . . , a s 0 ) is a weight factors vector, X(T 0 n) is the output time series, Y(T 0 n) is the FIR filter input-the hourly normalized Poisson data time series from MH according to (2), which begins from the values . .. Transfer function (TF) H(jωT 0 , a) for filter (2) is defined as follows: Here ω is the TF frequency parameter. For (3), a normalized fequence is introduced, w, ωT 0 = wπ, 0 ≤ w ≤ 1.0, and its discrete values are calculated: w l dw = 1.0/L 0 , w l = dw(l − 1), l = 1, . . . , L, L = L 0 + 1.
The frequency response (FR) H(w l , a) = |H(jw l , a)|, considering (3), is the following: for discrete normalized frequencies w l , l = 1, . . . , L according to (4). The cutoff frequency w c for FR (5) is found based on the equality |H(jw A , a)| 2 = 0.5. For a low-frequency FIR filter synthesis, the FR of the prototype filter is used, based on a Gaussian function H 0g (w, w c0 )

Synthesis Requiements
The problem of synthesis of the supposed FIR filter is solved based on the approximation of the FR function H 0g (w l , w c ) (6) in discrete points w l , l = 1, . . . , L with a FR function H g (w l , 0) according to (5). A functional S(H 0g , a, w c ) is formed: Obviously, the FR (5) represents a function which is polyharmonic in frequency w l . In case the prototype filter FR frequency derivative has discontinuities or is subject to strong alternations, e.g., if FR is a trapezoidal function, then the FR of the synthesized FIR filter, obtained based on approximation, will contain fluctuations due to the so-called Gibbs effect. Elimination and reduction of these fluctuations are usually achieved by choosing a suitable smooth prototype filter FR function. The smoothness requirement is largely satisfied by the Gaussian function (6).It should be noted that the Gaussian function is naturally suitable for the FR of a low-pass filter, since its values (6) practically differ from zero only in the region of low frequencies.
The requirements listed in the Introduction lead to formalized requirements: a.
Ensuring that the gain at zero frequencies is equal to one: Ensuring non-negativity of coefficients: For the synthesis procedure, it is assumed to set a small value r 0 , based on the a priori known duration of fluctuations, and some small cutoff frequency value w c for a prototype filter. The quasi-Gaussian filter synthesis procedure, consisting of finding the optimal coefficients a • s , s = 0, 1, . . . , s 0 , taking into account the requirements a,b, Equations (8) and (9) the predefined r 0 and w c , is performed on the basis of the approximation problem, which reduces to the implementation of conditional minimization: For a given small dimension r 0 of the synthesized quasi-Gaussian filter and a given small cutoff frequency w c for a prototype filter, the value for cutoff frequency to be found for a quasi-Gaussian filter is w cg , and the filter FR for the frequencies w l is denoted as The minimization of (10) could be performed based on modified direct zero-order optimization methods, taking into account the restrictions (8) and (9). However, because the (7) functional is multi-extremal, traditional modified direct methods, for example, using the coordinate descent method, the Hook-Jeeves method, the random descent method, etc. [13] do not provide successful minimization. The listed methods, as a rule, lead to "getting stuck" with search procedures in local minima.

Quasi-Gaussian Filter Synthesis Procedure
We can synthesize the quasi-Gaussian filter based on the technology of stochastic global minimization of the (7) functional with the constraints (8) and (9) using the optimization algorithm for annealing simulation [14,15]. To implement it, we will use the simulannealbnd.mat software module from the Matlab Global Optimization Toolbox [16].
Let us form a parallelepiped of constraints A 0 of dimension r 0 with boundaries a r , r = 1, . . . , r 0 -a ∈ A 0 , A 0 = {a : (0 ≤ a r ≤ a r , r = 1, . . . , r 0 )} and a new-with respect to (7)-functional S(H 0g , a) with a penalty term taking into account the constraint equality (8). Let us implement the global minimization of S(H 0g , a) taking into account A 0 using [16].
Let us set the initial vectors for the first iteration a 1 (I) ∈ A 0 , uniformly distributed in A 0 , I-a single descent procedure , I = 1, 2, . . . , I 0 , I 0 -a total number of descent procedures. Let us assume that each descent procedure consists of m 0 -a total number of iterations, m-a single iteration, m = 1, 2, . . . , m 0 . During descent, we assume that the initial value of the vector of parameters for (m + 1)-st iteration is equal to the calculated optimal value for the vector of parameters for m-th iteration-a m+1 (I) = a • m (I) . In each iteration, we perform n 0 descent steps, n is a descent step number, n = 1, 2, . . . , n 0 . Next, we will calculate the sequence of the functional S(m 0 , I)= S(H 0g , a • ) values and the corresponding optimal vectors a • (m 0 , I), I = 1, 2, . . . , I 0 . For global minimization, we search for the optimal index I • corresponding to the minimum of the S(m 0 , I) functional, and the optimal vector a • using brute force: On Figure 1, the example plots of the minimized S(m, I) functionals are displayed, depending on iteration number m and the descent procedure number I. Functionals are shown starting with m = 2, since for m = 1 their values are very large. Here, m 0 = 20; as the iteration number increases, the values of the functionals decrease. During the optimization process, a movement is made in a r 0 -dimensional space from one local minimum to another. Let us consider an example of quasi-Gaussian FIR filter synthesis. Based on the analysis of hourly experimental MH data from [12], it was found that the durations of possible fluctuations of the mathematical expectation in them were, on average, ≈10 ÷ 20 h and more. The dimension value r 0 , that could possibly allow the recognition of such fluctuations in mathematical expectations, was equal to 8. For a prototype filter FR (6), the parameter w c0 was related to the assigned cutoff frequency w c based on (6) (0.5) 1/2 = exp(−(w c /w c0 ) 2 ), w c0 = w c /(0.5 · ln 2) 1/2 .
We assign the cutoff frequency w c = 0.1, find w c0 and define H 0g (w, w c )-the prototype filter FR. By defining L we set the number of discrete normalized frequencies w l of calculations of the functional (7) for 0 ≤ w l ≤ 1.0, let us assume that L = 100 in our calculations. The polyharmonic FR function |H(jw)| (5) is formed from components performing 1, 2, . . . , s 0 fluctuations in this interval. For the accepted values L and r 0 , one period of the polyharmonic component with the maximum frequency corresponding to the number s 0 in (5), accounted for ≈15 sampling points of normalized frequencies w l , l = 1, . . . , L, which fully provided a fairly accurate calculation of the functional (7) necessary for direct search.
Let us calculate the vector of factors a • , form the synthesized quasi-Gaussian filter FR H g (w, w cg , a • ) and define the cutoff frequency w cg = 0.175 .
For the comparison, let us synthesize a common FIR filter using the fir1.mat module [3]. For the dimension r 0 = 8 and the assigned cutoff frequency w c = 0.1 we find out the final cutoff frequency w c f = 0.275; let us denote the FR as H f (w, w c f ). On Figure 2, the FR plots for H 0g (w, w c ), H g (w, w cg , a • ), H f (w, w c f ) are displayed. It is seen that, in case of low r 0 , the quasi-Gaussian filter FR was characterized by a better approximation to the prototype filter FR than the one of the common FIR filter.
Note that the proposed FIR filter, with the same dimension as the common FIR filter, made it possible to provide a lower value of the cutoff frequency than the realized cutoff frequency for the common FIR filter. The calculated cutoff frequencies of resulting FRs for common FIR filters synthesized using frequency sampling method and Fourier transforms [3,4] insignificantly (by ≈5-7% ) differ from the cutoff frequency w c f = 0.275. This gives a reason to make a conclusion about the advantages of a quasi-Gaussian filter over standard FIR filters.

Testing on Model Hourly Data Using Statistical Modeling
Testing of the proposed method and quasi-Gaussian filter was carried out on model hourly normalized data using statistical modeling [17]. For this purpose, on the basis of the Matlab module exprnd.mat [18], exponentially distributed model random numbers τ i , i = 1, 2, . . . were generated, with their mean value τ M0 , and the evenly distributed random registration time intervals T 0k , k = 1, 2, . . . within the range T 0 min ≤ T 0k ≤ T 0 max . The number of Poisson model events N M (Tk) was counted on the registration time intervals T 0k . Finding N M (Tk) was carried out by solving the conditional maximization problems: where for T 0k , k = 1, 2, . . . k 0 the range bounds T 0 min = 57 s, T 0 max = 59.5 s were assigned (see the Introduction section). Initial model 1-minutesampled and normalized Poisson-distributed data were constructed according to (11) The modulation of the average number of Poisson events in order to model decreases (increases) in the mathematical expectation was carried out by specifying the mean value function τ M0 (Tk) on the intervals (T(k − 1), Tk) for k from (12). For this, the relative modulation function µ(Tk), k = 1, 2, . . . , k 0 was formed and the initial temporal index of the modulation decrease k a , the duration of the decrease dk a and the depth of the relative decrease dµ ;. The function µ(Tk) was represented by the relations µ(Tk) = 1 − dµ for k a ≤ k ≤ k a + dk a , µ(Tk) = 1 for 1 ≤ k < k a , k a + dk a + 1 ≤ k ≤ k 0 . For the calculation example, the average number of Poisson model events per minute was set N M0 = 25, normalized average N M0 = N M0 /T, modulated normalized mean N M0 (Tk) = N M0 µ(Tk) = N M0 µ(Tk)/T, k = 1, 2, . . . , k 0 and the parameter τ M0 (Tk) = 1/(N M0 (Tk) − 1) was calculated.
For modeling, we assumed k 0 = 6000, which corresponded to the model minute data produced during 4.166 days. For the modulation function, the values k a = 1920, dk a = 1440 and dµ = 0.02 were taken. Model hourly averaged data Y M (T 0 n) for (13) with n 0 = 100, n a1 = 32, n a2 = n 1 + dn a dn a = 24. Figure 3 shows an example of statistical modeling results: the jagged light gray line with index 1 displays the Y E (T 0 n) plot; the solid line with index 2 denotes the fragment of X EG (T 0 n) which is the result of filtering the model dataset using a quasi-Gaussian filter; for comparison, the dashed line with index 3 denotes the fragment X EF (T 0 n) which is the result of filtering the model dataset using the software module fir1.mat [3]. Model piecewise constant modulating function Y M0 (T 0 n) = N M (T 0 n), represented by a dotted line (index 4), m 0 + dm = Y M0 (T 0 n) = 0.4165 for 1 ≤ n < n a1 , n a2 ≤ n < n 0 , m 0 = Y M0 (T 0 n) = 0.4087 for n a1 ≤ n ≤ n a2 , where the value of dm = 0.833 × 10 −2 corresponded to the predefined 2% decrease. The plots show that the result of the quasi-Gaussian filter application (line 2) is a better approximation to the model piecewise constant modulation (line 4) than the result of a common FIR filtering (line 3).  The calculation of approximate estimates of filtering errors for the quasi-Gaussian filter and fir1-filter was performed by calculating the root-mean-square (RMS) errors according to the following formulas for datasets Y M0 (T 0 n), Y MG (T 0 n), Y MF (T 0 n): Results of a large number of tests performed for (14) showed that the σ MG error values for X MG (T 0 i) regarding Y MG (T 0 n) are, on average, 15-30% less than the corresponding σ MF error values for X MF (T 0 i). An overview of model X MG (T 0 i) and X MF (T 0 i) (Figure 3) made it possible to ensure that the minimum duration of the interval, within which recognition for the decrease dµ = 0.02 can be performed, is 12-24 h.
The proposed method and the quasi-Gaussian filter provided more noise reduction than a common FIR filter. Consideration of the results of statistical modeling made it possible to draw a conclusion about the efficiency of the quasi-Gaussian filtering method.

Testing the Method and the Quasi-Gaussian Filter on Experimental Normalized Poisson Data from the URAGAN Hodoscope
Testing in this section consisted of determining the performance and capabilities of the proposed method and the quasi-Gaussian filter for recognizing small in duration and magnitude decreases in time intervals for the experimental hourly normalized Poisson data registered by the URAGAN hodoscope, taken from [12].
For analysis, a time interval was selected from 09/02/2017, 20:00 UTC to 09/18/2017, 15:00 UTC, with a total duration of 15.6 days. During this interval, the heliosphere was turbulent due to strong solar coronal mass ejections (CMEs) The CMEs that occurred on that period, caused intense geomagnetic storms that were discussed, for example, in [19,20]. The emerging CMEs caused modulations of muon fluxes recorded in MH and led to lower mathematical expectations (including the ones due to Forbush decreases) in Poisson MH data.
MH data were the matrix series of distribution functions of the intensities of muon fluxes Y E (i, j, T 0 n), defined in a rectangular region i = 1, . . . ,N 1 , j = 1, . . . ,N 2 , N 1 = 90, N 2 = 76, n = 1, 2, . . .. Solid angles correspond to azimuth and zenith indices i, j, ϕ i = ∆ϕ(i − 1), ϑ j = ∆ϑ(j − 1) , ∆ϕ = 1 • , ∆ϑ = 4 • in which the registered particles were counted. MH data Y E (j 0 , i 0 , T 0 n) were a time series with indices j 0 , i 0 ; the considered interval was determined for n E min ≤ n ≤ n E max , n E min = 5900, n E max = 6275 (counting hours for [12] began from the first hour of 2017). Figure 5 shows the results of quasi-Gaussian filtering and interval recognition with reductions in mathematical expectation. The original data Y E (T 0 n) for j 0 = 30, i 0 = 31 were denoted by light gray jagged lines (index 1). Fluctuations in data with a period of ≈24 h and an amplitude of ≈0.0037-0.0040 are due to the daily rotation of the MH with the Earth. Line with index 2 depicts the data X EG (T 0 n) filtered based on quasi-Gaussian filter. The recognized intervals of intensity decrease, intensity recovery and intensity mathematical expectation decrease were denoted by a piecewise linear spline-like dashed line X ES (T 0 n) (index 3). Analysis of intervals 5969-6043, 6127-6189 based on X ES (T 0 n) leads to a conclusion that the mathematical expectation values of decreases on them were ∆m 1 = 0.01, ∆m 2 = 0.005 for the relative decrease rates dµ 1 = 0.027, dµ 2 = 0.020. For Y E (T 0 n) and X EG (T 0 n), the mathematical expectations on these intervals were m • E1 = 0.3505, m • E1 = 0.3510, and m • EG1 = 0.490, m • EG2 = 0.3520, respectively, on average. The errors of the mathematical expectations estimates were ∆m • = 0.0010 − −0.0015, which is 10-30% from the mathematical expectation values obtained, and this led to successful recognition of decreases with the relative decrease rates of 0.02-0.03.
Testing on experimental MH data made it possible to draw a conclusion about the efficiency of the quasi-Gaussian filtering method and its satisfactory capabilities for recognizing small fluctuations of the mathematical expectations.

Discussion
The comparison between the model data filtering result obtained using the proposed filter and the one obtained using the fir1 (plots on Figure 3) shows that the resulting time series are close to each other; however, the X MG (T 0 i) seems to be closer to the initial model. The main quantitative result of testing the method and the quasi-Gaussian filter on model normalized Poisson datasets included the calculations for (14) for a set of realizations/ The resulting errors σ MG for X MG (T 0 i), on average, by 15-30% less errors σ MF for X MF (T 0 i). This means that the proposed filtering method provided better filtering (noise reduction) than the standard FIR filter. Consideration of the results of statistical modeling made it possible to draw a conclusion about the efficiency of the method and the quasi-Gaussian filter.
Further tests of the new method on model data, aimed at estimating the mathematical expectation and its RMS errors with respect to the durations and magnitudes of model decreases, showed the method capabilities in disturbance recognitions. It can be seen on Figure 4 that the coefficients ε • gm turned out to be less than the values of the coefficients ε • f m , on average, by about 10-30%. The nature of the plots of coefficients ε • gσ and ε • f σ for the RMS for the same parameters dn a , dµ is almost the same-the coefficients ε • gσ are less than the values of the coefficients ε • f σ , on average, also by ≈10-30%. From the point of view of recognizing decreases in duration and magnitude, the use of a quasi-Gaussian filter is more preferable than a common FIR filter.
Finally, tests made on real experimental datasets from a muon hodoscope display the method application to data processing and recognition of intervals of decreasing and recovering muon flux intensity. Due to the noise reduction in X EG (T 0 n), it became possible to clearly see the intervals of quiet data ( Figure 5), intervals with decreases and recoveries and intervals with declines in mathematical expectation; all these recognized intervals were denoted by a line X ES (T 0 n) (index 3 on Figure 5). On the intervals with the boundary points 5900-5954, 6057-6121, 6197-6276 there were quiet data, on the time intervals 5969-6043, 6127-6189 a decrease in mathematical expectation was observed, the time intervals 5955-5970, 6044-6056, 6122-6126, 6190-6196 corresponded to data with decreases and recoveries. On the intervals 5969-6043, 6127-6189, it is quite possible to recognize relative reductions in mathematical expectation. The errors of the mathematical expectations estimates were ∆m • = 0.0010-0.0015, which is 10-30% from the mathematical expectation values obtained, and this led to successful recognition of decreases with the relative decrease rates of 0.02-0.03 and an average duration of mathrm approx 10 h.
Testing the proposed method and quasi-Gaussian filter for data variants with indices j 0 = 31, i 0 = 30, allowed to obtain results that are almost similar to those depicted on Figure 5); the errors in the estimation of the boundary points of the sections during recognition with depressions amounted to δn ≈ 2-5 h. Thus, testing on experimental MH data allowed us to make a conclusion about the efficiency of the method and the quasi-Gaussian filter and their satisfactory capabilities for recognizing mathematical expectation small in duration and magnitude.

Conclusions
The proposed filtering method for time series of normalized Poisson-distributed data, which was based on the developed digital low-pass quasi-Gaussian filter with a finite impulse response, a gain equal to one at low frequencies and non-negative weighting coefficients, turned out to be efficient; the FR of the low-frequency quasi-Gaussian filter of small dimension was characterized by a better approximation to the prototype filter FR than the FR of common FIR filters.
Testing the filtering method based on the quasi-Gaussian filter for the problems of recognizing small in duration and magnitude fluctuation decreases (increases) in mathematical expectations using statistical modeling and statistical tests have confirmed its effectiveness: • The proposed method provided a decrease in errors in the filtered time series in comparison with the error values for standard FIR filters, by ≈15-30%; the method made it possible to recognize the mathematical expectation fluctuations with a relative decrease of 0.02 and duration of ≈12-24 h; • The proposed method and the developed quasi-Gaussian filter provided relative error coefficients for mathematical expectation and root mean square values that appeared to be ≈10-30% less than the error coefficients for common FIR filters. Testing the method and the low-frequency quasi-Gaussian filter on experimental Poisson data made it possible to draw a conclusion about its satisfactory capabilities for recognizing decreases with relative decrease coefficients ≈0.020-0.030.
The proposed method of noise reduction and a quasi-Gaussian filter have favorable prospects of using radiation particle counters for digital information processing in problems of experimental physics and measuring technology. Data Availability Statement: Data sharing is not applicable to this article.