Next Article in Journal
Exploring Hybrid-Multimodal Routing to Improve User Experience in Urban Trips
Next Article in Special Issue
The August 2018 Geomagnetic Storm Observed by the High-Energy Particle Detector on Board the CSES-01 Satellite
Previous Article in Journal
Exploitation of Antimicrobial Nanoparticles and Their Applications in Biomedical Engineering
Previous Article in Special Issue
Trapped Proton Fluxes Estimation Inside the South Atlantic Anomaly Using the NASA AE9/AP9/SPM Radiation Models along the China Seismo-Electromagnetic Satellite Orbit
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Low-Pass Filtering Method for Poisson Data Time Series

1
Geophysical Center of the Russian Academy of Sciences (GC RAS), 119296 Moscow, Russia
2
Schmidt Institute of Physics of the Earth of the Russian Academy of Sciences (IPE RAS), 123242 Moscow, Russia
3
Scientific & Educational Centre NEVOD, National Research Nuclear University MEPhI (NRNU MEPhI), 115409 Moscow, Russia
*
Authors to whom correspondence should be addressed.
Submission received: 23 April 2021 / Revised: 11 May 2021 / Accepted: 12 May 2021 / Published: 15 May 2021
(This article belongs to the Special Issue Research on Cosmic Rays and Their Impact on Human Activities)

Abstract

:
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 filtering method for normalized Poisson-distributed data time series is proposed. A digital quasi-Gaussian filter is designed, with a finite impulse response and non-negative weights. The quasi-Gaussian filter synthesis is implemented using the technology of stochastic global minimization and modification of the annealing simulation algorithm. The results of testing the filtering method and the quasi-Gaussian filter on model and experimental normalized Poisson data from the URAGAN muon hodoscope, that have confirmed their effectiveness, are presented.

1. 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 ( T k ) —the quantities of Poisson-distributed events recorded in a given solid angle at time intervals ( T ( k 1 ) , T ( k 1 ) + T 0 k ) , k = 1 , 2 , , k 0 , where T = 1 minute. Due to the features of the MH design, registration intervals T 0 k are random with a uniform distribution law in the range T 0 min T 0 k T 0 max < T .
From each secondary MH transducer, the 1-minute-sampled normalized Poisson data Y ( T k ) 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:
Y ( T k ) = N ( T k ) / T 0 k , Y ( T 0 n ) = 1 60 k = 1 + 60 ( n 1 ) k = 60 n Y ( T k ) , n = 1 , 2 , , T 0 = T · 60 .
Data resulting from (1) are produced for the whole set of solid angles; next, they are placed into time series of matrix MH data in the database [12].

2. Method

2.1. Quasi–Gaussian Digital Low-Pass Filter

2.1.1. Statement of the Problem

One-dimensional FIR filter synthesized here is built according to the following difference equation:
X ( T 0 n ) = s = 0 s 0 a s Y ( T 0 ( n s ) ) , n = 1 , 2 , ,
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 Y ( T 0 ( 1 s 0 ) ) , Y ( T 0 ( 1 s 0 + 1 ) ) , Y ( T 0 ( 1 s 0 + 2 ) ) , . Transfer function (TF) H ( j ω T 0 , a ) for filter (2) is defined as follows:
H ( j ω T 0 , a ) = s = 0 s 0 a s e j 2 π ω T 0 s .
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
d w = 1.0 / L 0 , w l = d w ( l 1 ) , l = 1 , , L , L = L 0 + 1 .
The frequency response (FR) H ( w l , a ) = H ( j w l , a ) , considering (3), is the following:
H ( w l , a ) 2 = H 1 2 ( w l , a ) + H 2 2 ( w l , a ) , H 1 ( w l , a ) = s = 0 s 0 a s cos ( 2 π w l s ) , H 2 ( w l , a ) = s = 0 s 0 a s sin ( 2 π w l s )
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 ( j w 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 0 g ( w , w c 0 )
H 0 g ( w , w c 0 ) = exp ( ( w / w c 0 ) 2 ) .

2.1.2. Synthesis Requiements

The problem of synthesis of the supposed FIR filter is solved based on the approximation of the FR function H 0 g ( 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 0 g , a , w c ) is formed:
S ( H 0 g , a , w c ) = l = 1 L [ ( s = 0 s 0 a s C s ( w l ) ) 2 + ( s = 0 s 0 a s S s ( w l ) ) 2 H 0 g 2 ( w l , w c ) ] 2 .
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:
  • Ensuring that the gain at zero frequencies is equal to one:
    1 = H ( 0 , a ) = s = 0 s 0 a s , a A 1 , A 1 = { a : ( 1 = s = 0 s 0 a s ) } ;
  • Ensuring non-negativity of coefficients:
    a A 0 , A 0 = { a : ( 0 a s , s = 0 , 1 , , s 0 ) } ;
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:
a ( w c ) = arg { min a A 0 , a A 1 S ( H 0 g , a , w c ) } .
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 c g , and the filter FR for the frequencies w l is denoted as H g ( w l , w c g , a ) , l = 1 , , L .
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.

2.2. 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 0 g , a ) with a penalty term taking into account the constraint equality (8). Let us implement the global minimization of S ¯ ( H 0 g , 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 0 g , 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:
I = arg ( min 1 I I 0 S ( m 0 , I ) } , a = a ( m 0 , I ) .
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 c 0 was related to the assigned cutoff frequency w c based on (6)
( 0.5 ) 1 / 2 = exp ( ( w c / w c 0 ) 2 ) , w c 0 = w c / ( 0.5 · ln 2 ) 1 / 2 .
We assign the cutoff frequency w c = 0.1 , find w c 0 and define H 0 g ( 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 ( j w ) (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 c g , a ) and define the cutoff frequency w c g = 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 0 g ( w , w c ) , H g ( w , w c g , 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.

3. Results

3.1. Testing the Method and the Quasi–Gaussian Filter on Model Normalized Poisson Data

3.1.1. 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 τ M 0 , and the evenly distributed random registration time intervals T 0 k , k = 1 , 2 , within the range T 0 m i n T 0 k T 0 m a x . The number of Poisson model events N M ( T k ) was counted on the registration time intervals T 0 k . Finding N M ( T k ) was carried out by solving the conditional maximization problems:
N M ( T k ) = arg { max N M } i > 0 ,
providing that T 0 k i = 1 N M τ , where for T 0 k , k = 1 , 2 , k 0 the range bounds T 0 m i n = 57 s, T 0 m a x = 59.5 s were assigned (see the Introduction section). Initial model 1-minute-sampled and normalized Poisson-distributed data were constructed according to (11) and the calculation of relations N ¯ M ( T k ) , similar to (1):
N ¯ M ( T k ) = N M ( T k ) / T 0 k , k = 1 , 2 , k 0 .
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 τ M 0 ( T k ) on the intervals ( T ( k 1 ) , T k ) for k from (12). For this, the relative modulation function μ ( T k ) , k = 1 , 2 , , k 0 was formed and the initial temporal index of the modulation decrease k a , the duration of the decrease d k a and the depth of the relative decrease d μ ;. The function μ ( T k ) was represented by the relations μ ( T k ) = 1 d μ for k a k k a + d k a , μ ( T k ) = 1 for 1 k < k a , k a + d k a + 1 k k 0 .
For the calculation example, the average number of Poisson model events per minute was set N M 0 = 25 , normalized average N ¯ M 0 = N M 0 / T , modulated normalized mean N ¯ M 0 ( T k ) = N ¯ M 0 μ ( T k ) = N M 0 μ ( T k ) / T , k = 1 , 2 , , k 0 and the parameter τ M 0 ( T k ) = 1 / ( N ¯ M 0 ( T k ) 1 ) was calculated.
Based on [18], random exponentially distributed numbers with τ M 0 ( T k ) and random evenly distributed values with T 0 m i n = 57 s , T 0 m a x = 59 , 5 s were generated, with the use of which by (11), model Poisson data N M ( T k ) and by (12)—normalized Poisson data N ¯ M ( T k ) were calculated. Further, similarly to (1), a time series of averaged model hourly normalized Poisson data was formed:
Y M ( T 0 n ) = 1 60 k = 1 + 60 ( n 1 ) k = 60 n N ¯ M ( T k ) , n = 1 , 2 , , n 0 , n 0 = e n t ( k 0 / 60 ) .
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 , d k a = 1440 and d μ = 0.02 were taken. Model hourly averaged data Y M ( T 0 n ) for (13) with n 0 = 100 , n a 1 = 32 , n a 2 = n 1 + d n a d n 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 E G ( 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 E F ( 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 M 0 ( T 0 n ) = N ¯ M ( T 0 n ) , represented by a dotted line (index 4), m 0 + d m = Y M 0 ( T 0 n ) = 0.4165 for 1 n < n a 1 , n a 2 n < n 0 , m 0 = Y M 0 ( T 0 n ) = 0.4087 for n a 1 n n a 2 , where the value of d m = 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 M 0 ( T 0 n ) , Y M G ( T 0 n ) , Y M F ( T 0 n ) :
σ M G 2 = 1 n 0 n = 1 n 0 ( Y M 0 ( T 0 n ) Y M G ( T 0 n ) ) 2 , σ M F 2 = 1 n 0 n = 1 n 0 ( Y M 0 ( T 0 n ) Y M F ( T 0 n ) ) 2 .
Results of a large number of tests performed for (14) showed that the σ M G error values for X M G ( T 0 i ) regarding Y M G ( T 0 n ) are, on average, 15–30% less than the corresponding σ M F error values for X M F ( T 0 i ) . An overview of model X M G ( T 0 i ) and X M F ( 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.

3.1.2. Estimation of Mathematical Expectation and Its Root Mean Square Errors

Testing of the method and quasi-Gaussian filter for estimating the mathematical expectation and the RMS of its errors depending on d n a —the duration of decreases and d μ —the relative decrease value were carried out using statistical tests [17]. Random datasets Y M ( s , T 0 i ) , X M G ( s , T 0 i ) , X M F ( s , T 0 i ) , s = 1 , 2 , , M , where s is the number of the dataset, M is the total quantity of datasets. The estimates of mathematical expectation m g ( d n a , d μ ) and RMS values σ g ( d n a , d μ ) for X M G ( s , T 0 i ) for a set of values d n a and d μ
m g ( s , d n a , d μ ) = 1 n a n = n a 1 n a 1 + d n a X M G ( s , T 0 n ) , m g ( d n a , d μ ) = 1 M s = 1 M m g ( s , d n a , d μ ) , σ g ( s , d n a , d μ ) = 1 n a 1 n = n a n a + d n a ( X M G ( s , T 0 n ) m g ( s , d n a , d μ ) ) 2 , σ g ( d n a , d μ ) = 1 M s = 1 M σ g ( s , d n a , d μ ) .
The coefficients of relative errors ε g m ( d n a , d μ ) , ε g σ ( d n a , d μ ) of the quasi-Gaussian filter as ratios of errors m g ( d n a , d μ ) m 0 and RMS σ g ( s , d n a , d μ ) to the values of d m reductions are the following:
ε g m ( d n a , d μ ) = ( m g ( d n a , d μ ) m 0 ) / d m , ε g σ ( d n a , d μ ) = ( σ g ( d n a , d μ ) ) / d m
The coefficients ε g m , ε g σ , calculated for d n a , d μ , characterized the recognition capabilities of quasi-Gaussian filtering model decreases. Similarly, using (15) and (16) m f ( d n a , d μ ) and σ f ( d n a , d μ ) for X M F ( s , T 0 n ) and the coefficients ε f m ( d n a , d μ ) , ε f σ ( d n a , d μ ) . On Figure 4, the results of statistical tests are displayed, where M = 500 . The ε g m ( d n a , d μ ) coefficients plots are the solid lines with indices 1, 2, and the ε f m ( d n a , d μ ) plots are the dashed lines with indices 3, 4. The coefficients ε g m , ε f m are given depending on the duration with the values d n a = 12 , 24 , 48 , 72 h and relative decreases in d μ , taking the values of 0.01 (indices 1, 3) and 0.03 (indices 2, 4).
The effect of quasi-Gaussian filtering was determined based on the calculation of δ ε f g , m —the rates of errors with respect to the mathematical expectations:
δ ε f g , m ( d n a , d μ ) = ( ε f m ( d n a , d μ ) ε g m ( d n a , d μ ) ) / ε g m ( d n a , d μ )
The results of the δ ε f g , m calculations according to (17) for some d μ and d n a values are:
  • δ ε f g , m = 0.115 (11.5%) for d n a = 24 and d μ = 0.01 ;
  • δ ε f g , m = 0.196 (19.6%) for d n a = 24 and d μ = 0.03 .
Analysis of the error values showed that the ε g m rate values appeared to be about 10–30% lower than the ε f m values. The nature of the dependencies of the estimates of the error coefficients for the ε g σ and ε f σ root mean square values for the same d n a and d μ parameters is almost the same: the ε g σ are also ≈10–30% lower than the ε f σ . This means that, for the recognition of decreases small in duration and magnitude, the use of a quasi-Gaussian filter is more preferable than the use of a common FIR filter.

3.2. 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 E G ( 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 E S ( T 0 n ) (index 3). Analysis of intervals 5969–6043, 6127–6189 based on X E S ( 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 E G ( T 0 n ) , the mathematical expectations on these intervals were m E 1 = 0.3505 , m E 1 = 0.3510 , and m E G 1 = 0.490 , m E G 2 = 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.

4. 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 M G ( 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 σ M G for X M G ( T 0 i ) , on average, by 15–30% less errors σ M F for X M F ( 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 ε g m 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 d n 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 E G ( 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 E S ( 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 m a t h r m a p p r o x 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.

5. 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.

Author Contributions

Conceptualization, V.G., A.S. and I.Y.; methodology, V.G.; software, V.G. and V.C.; validation, M.D., A.G. and A.S.; formal analysis, V.C.; investigation, A.S.; resources, A.G.; data curation, A.D., A.K. and N.O.; writing—original draft preparation, V.G.; writing—review and editing, R.S. and M.D.; visualization, V.G.; supervision, I.Y.; project administration, A.S.; funding acquisition, A.G. All authors have read and agreed to the published version of the manuscript.

Funding

This work was funded by the Russian Science Foundation (project No. 17-17-01215).

Data Availability Statement

Data sharing is not applicable to this article.

Acknowledgments

The results of experiments presented in this research rely on data collected by the Scientific & Educational Centre NEVOD, National Research Nuclear University MEPhI. We acknowledge URAGAN muon hodoscope data provided by the NEVOD institution. This work employed facilities and data provided by the Shared Research Facility “Analytical Geomagnetic Data Center” of the Geophysical Center of RAS (http://ckp.gcras.ru/) (accessed on 19 April 2021). We would like to thank two anonymous reviewers whose valuable comments helped to improve the manuscript and properly demonstrate the results of our research.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Lloyd, E.; Ledermann, W. (Eds.) Handbook of Applicable Mathematics, Russian ed.; Finansy I Statistika: Moscow, Russia, 1989; Volume 1, pp. 45–83. (In Russian) [Google Scholar]
  2. Taylor, F.J. Digital Filters: Principles and Applications with MATLAB; J. Wiley I& Sons: New York, NY, USA, 2011; pp. 53–70. [Google Scholar]
  3. Filter Design Matlab Toolbox. Available online: http://matlab.exponenta.ru (accessed on 7 March 2021).
  4. Sergienko, A.B. Digital Signal Processing: Textbook, 3rd ed.; BHV-Peterburg: St. Petersburg, Russia, 2011; pp. 371–440. (In Russian) [Google Scholar]
  5. Getmanov, V.G. Digital Signal Processing, 2nd ed.; National Research Nuclear University MEPhI: Moscow, Russia, 2010; pp. 210–225. (In Russian) [Google Scholar]
  6. Goncharov, G.A.; Zubatkin, O.Y.; Lopatin, P.A. Calculation of a filter with a frequency response close to a Gaussian curve. Commun. Equip. Dig. 1978, 5, 43–49. (In Russian) [Google Scholar]
  7. Bugrov, V.N.; Bessonova, E.V. Numerical design of Gaussian digital filters. Electron. Des. Technol. 2012, 3, 40–42. (In Russian) [Google Scholar]
  8. Rockenbach, M.; Dal Lago, A.; Schuch, N.J.; Munakata, K.; Kuwabara, T.; Oliveira, A.G.; Echer, E.; Braga, C.R.; Mendonça, R.R.S.; Kato, C.; et al. Global Muon Detector Network Used for Space Weather Applications. Space Sci. Rev. 2014, 182, 1–18. [Google Scholar] [CrossRef]
  9. Grupen, C.; Shwartz, B. Particle Detectors, 2nd ed.; Cambridge University Press: New York, NY, USA, 2008; pp. 56–84. [Google Scholar]
  10. Yashin, I.I.; Astapov, I.I.; Barbashina, N.S.; Borog, V.V.; Chernov, D.V.; Dmitrieva, A.N.; Kokoulin, R.P.; Kompaniets, K.G.; Mishutina, Y.N.; Petrukhin, A.A.; et al. Real-time data of muon hodoscope URAGAN. Adv. Space Res. 2015, 56, 2693–2705. [Google Scholar] [CrossRef]
  11. Barbashina, N.S.; Kokoulin, R.P.; Kompaniets, K.G.; Mannocchi, A.; Petrukhin, A.A.; Timashkov, D.A.; Saavedra, O.; Trinchero, G.; Chernov, D.V.; Shutenko, V.V.; et al. The URAGAN wide-aperture large-area muon hodoscope. Instrum. Exp. Technol. 2008, 51, 180–186. [Google Scholar] [CrossRef]
  12. NEVOD COMPLEX. National Research Nuclear University MEPhI. Available online: http://www.nevod.mephi.ru (accessed on 6 March 2021).
  13. Himmelblau, D.M. Applied Nonlinear Programming, Russian ed.; Mir: Moscow, Russia, 1975; pp. 157–193. (In Russian) [Google Scholar]
  14. Panteleev, A.V.; Skavinskaya, D.V. Metaheuristic Algorithms for Global Optimization; University Book: Moscow, Russia, 2019; pp. 5–29. (In Russian) [Google Scholar]
  15. Ingber, L.; Oliveira, E.H.; Petraglia, A.L.; Petraglia, M.R.; Machado, M.A.S. Stochastic Global Optimization and its Applications with Fuzzy Adaptive Simulated Annealing; Springer: Berlin/Heidelberg, Germany, 2012; pp. 33–62. [Google Scholar]
  16. Global Optimization Matlab Toolbox. Available online: http://matlab.exponenta.ru (accessed on 4 March 2021).
  17. Mikhailov, G.A.; Voytishek, A.V. Numerical Statistical Modeling. Monte Carlo Method; Yurayt Publishing House: Moscow, Russia, 2018; pp. 126–174. (In Russian) [Google Scholar]
  18. Statistic Matlab Toolbox. Available online: http://matlab.exponenta.ru (accessed on 5 March 2021).
  19. Sidorov, R.; Soloviev, A.; Gvishiani, A.; Viktor Getmanov, V.; Mandea, M.; Petrukhin, A.; Yashin, I.; Obraztsov, A. A combined analysis of geomagnetic data and cosmic ray secondaries for the September 2017 space weather event studies. Russ. J. Earth Sci. 2019, 19, ES4001. [Google Scholar] [CrossRef] [Green Version]
  20. Oshchenko, A.A.; Sidorov, R.V.; Soloviev, A.A.; Solovieva, E.N. Overview of anomality measure application for estimating geomagnetic activity. Geophys. Res. 2020, 21, 51–69. (In Russian) [Google Scholar]
Figure 1. Plots of descent procedures—minimization of functionals S ( m , I ) , I = 1 , 2 , , I 0 , m = 1 , 2 , , m 0 .
Figure 1. Plots of descent procedures—minimization of functionals S ( m , I ) , I = 1 , 2 , , I 0 , m = 1 , 2 , , m 0 .
Applsci 11 04524 g001
Figure 2. FR plots: H 0 g ( w , w c ) (line 1), H g ( w , w c g , a ) (line 2), H f ( w , w c f ) (line 3).
Figure 2. FR plots: H 0 g ( w , w c ) (line 1), H g ( w , w c g , a ) (line 2), H f ( w , w c f ) (line 3).
Applsci 11 04524 g002
Figure 3. Fragments of model datasets Y M ( T 0 n ) (line 1), filtering results X M G ( T 0 n ) (line 2), X M F ( T 0 n ) (line 3) and model modulating function Y M 0 ( T 0 n ) (line 4).
Figure 3. Fragments of model datasets Y M ( T 0 n ) (line 1), filtering results X M G ( T 0 n ) (line 2), X M F ( T 0 n ) (line 3) and model modulating function Y M 0 ( T 0 n ) (line 4).
Applsci 11 04524 g003
Figure 4. Results of calculating the coefficients of relative errors ε g m , ε f m .
Figure 4. Results of calculating the coefficients of relative errors ε g m , ε f m .
Applsci 11 04524 g004
Figure 5. Results of quasi-Gaussian filtering and identification of regions with Forbush decreases: original data (line 1), filtered data (line 2), recognized intervals of various muon flux intensity (line 3).
Figure 5. Results of quasi-Gaussian filtering and identification of regions with Forbush decreases: original data (line 1), filtered data (line 2), recognized intervals of various muon flux intensity (line 3).
Applsci 11 04524 g005
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Getmanov, V.; Chinkin, V.; Sidorov, R.; Gvishiani, A.; Dobrovolsky, M.; Soloviev, A.; Dmitrieva, A.; Kovylyaeva, A.; Osetrova, N.; Yashin, I. Low-Pass Filtering Method for Poisson Data Time Series. Appl. Sci. 2021, 11, 4524. https://0-doi-org.brum.beds.ac.uk/10.3390/app11104524

AMA Style

Getmanov V, Chinkin V, Sidorov R, Gvishiani A, Dobrovolsky M, Soloviev A, Dmitrieva A, Kovylyaeva A, Osetrova N, Yashin I. Low-Pass Filtering Method for Poisson Data Time Series. Applied Sciences. 2021; 11(10):4524. https://0-doi-org.brum.beds.ac.uk/10.3390/app11104524

Chicago/Turabian Style

Getmanov, Victor, Vladislav Chinkin, Roman Sidorov, Alexei Gvishiani, Mikhail Dobrovolsky, Anatoly Soloviev, Anna Dmitrieva, Anna Kovylyaeva, Nataliya Osetrova, and Igor Yashin. 2021. "Low-Pass Filtering Method for Poisson Data Time Series" Applied Sciences 11, no. 10: 4524. https://0-doi-org.brum.beds.ac.uk/10.3390/app11104524

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop