Reconstruction of digital filter parameters when changing the data arrival period

Introduction

In digital control of technology, radio communication and automation, digital filters (DF) are widely used, which provide processing of discrete samples of a continuous signal x(t) at fixed times t=k*T0, where Т0 (second) is the certain step of discreteness in time, k=0,1,2.. is the integer variable that defines the dimensionless discrete time. The implementation of such devices is usually carried out on circuits containing registers and adders with certain coefficients; implementation in the form of a digital microprocessor device is possible (Serrezuela et al., 2017; Chen et al., 2020). A certain filter functional is defined by a specific mathematical model, difference equation, by a corresponding z-transfer function, or by a discrete impulse transient function of the filter (Shamrikov, 1985; Oppenheim & Schafer, R. 2018; Antonio, 1983). It should be noted that these models determine the operation of the filter at a specific value of the step T0=1/f0, where f0 (hertz) is the frequency of periodic sample arrivals of x[k] continuous signal x(t) (Gadzikovsky, 2007; Adamou-Mitiche & Mitiche, 2019). The time intervals of signal arrival in many technical problems can change significantly during operation, for example, as it is described in (Psiaki et al., 2014; Daneshmand et al., 2015). In particular, this may be due to a change in the nature of transmission over the radio path (Velikanova & Voroshilin, 2012; Xu et al., 2016) and also due to the conditions for receiving and transmitting information from continuous objects. So, in works (Sokolov et al., 2018; Kalmykov et al., 2020) the changes associated with the satellite measurement system are considered. Without dwelling on the reasons for these phenomena, we will assume that the discreteness period T0=1/f0 during system operation can change to the value TN=N/f0, where N is some number characterizing the multiplicity of changes in the period relative to the calculated one (Koshita et al., 2017). If the frequency changes during operation f0 of data arrival on frequency fN=f0/N the properties of the filter will change, and with a significant change in frequency, such filtering can lead to unsatisfactory results, if the device parameters are not rebuilt. Particularly critical change of f0 is for digital automatic systems (Shamrikov, 1985), where DF operates in a closed control loop, and a change in their properties can lead to a decrease in stability margins, and, possibly, to a loss of stability of the digital system. We show it on the example of a linear digital filter of the n-th order. Let the processing of a continuous signal x(t) by a discrete sample x[k] is performed by a linear DF determined by z-transfer function of the following form:

where 𝑧 = 𝑒𝑠∙𝑇𝑜, s is the Laplace transform parameter, 𝐴(𝑧−1), 𝐵(𝑧−1) are the polynomials from z−1, determined by a specific set of coefficients 𝑎𝑖, 𝑏𝑖, i=0..n. DF of the described type, operating at a frequency f0 =1/T0, will be called the reference DF. The corresponding difference equation written with respect to the input y[k] and output x[k] of DF for the periods of time 𝑡 = 𝑘 ∙ 𝑇0 will take the following form:

This difference equation defines one of the ways to implement DF in the form of an algorithm for calculating the output coordinate x[k] by past output x[k-i] and input y[k-i] values detainees, i.e. saved in memory at the previous steps of calculation. For fixed frequency f0 DF with z-transfer function (1) performs a certain transformation of the spectrum of a continuous signal, for example, it filters out low-frequency or high-frequency signal components or components of a certain frequency, corrects the phase characteristic, etc. (Antonio, 1983; Oppenheim & Schafer, 2018). In this case, Fourier transform of the signal at the output of DF is determined by the following expression:

It should be noted that due to the specific properties of the transcendental function 𝑒𝑗∙𝜔∙𝑇𝑜 , spectrum 𝑋( 𝑒𝑗∙𝜔∙𝑇𝑜) of discrete signal x[k] is periodic with a period 𝜔𝑜(𝑟adians/second) = 2 ∗ 𝜋 ∗ 𝑓0 (Shamrikov, 1985). The samples are received not with the calculated interval T0, but with a changed one: TN=T0*N, where the parameter N determines the multiplicity of the change in the quantization period TN relative to the calculated period T0. An important point is that the new quantization period T0*N should also ensure the transmission of the properties of a continuous signal by discrete samples, which, according to the Nyquist– Shannon theorem, corresponds to the relation 𝜔𝑚 < 𝜋 𝑇0∗𝑁, where 𝜔𝑚(𝑟𝑎𝑑𝑖𝑎𝑛𝑠/𝑠𝑒𝑐𝑜𝑛𝑑) is the maximum frequency present in the spectrum of a continuous signal. When changing the quantization period TN=T0*N, the frequency response (FR) of DF will also change, which is determined by the following formula:

This FR has properties that differ from the properties of the reference filter: it is scaled along the frequency axis, and most importantly, it has a repetition period 2𝜋 ∙ 𝑓0⁄𝑁, which is N times less (if N>1), or more than the original one. For this reason, to choose a "similar filter", i.e. some DF DN (z, T0⋅ 𝑁) of the same order, which would completely repeat the properties of the reference DF with a change in frequency f0 to fN=f0/𝑁, is impossible. With regard to the simplest DF of 1-2 orders of magnitude, a possible approximate approach for generating a conversion is given in (Belonogov, 2014; Kaplun et al., 2020). In this paper, we pose the problem of calculating the algorithms for tuning the parameters of such a "similar" DF operating at the frequency fN, which would have closeness of its FR to the characteristics of the reference filter 𝐷0(𝑒𝑗∙𝜔∙𝑇𝑜) in the baseband of a digital system: 0 < 𝜔 < 𝜋 ∙ 𝑓0/𝑁. The procedure should ensure the constancy of the dynamic characteristics of DF, while allowing the use of algorithms in the operational restructuring of DF parameters when the period T0 changes. On the basis of studies carried out in the frequency domain using the methods of w-transformation, an algorithm for the operational restructuring of the numerical parameters of the filter based on information about the time intervals of information receipt has been developed. A method for calculating the restructuring of parameters based on the preliminary formation of a special recalculation matrix is proposed. At the stage of the filter operation in real time, this matrix implements the operational recalculation of the filter parameters by a linear relation. The proposed approaches make it possible to ensure the constancy of the frequency properties of the digital filter with high accuracy with a significant change in the periods of information receipt. At the same time, the implementation of the proposed approaches does not require significant computing resources and can be carried out on 8-16 bit microprocessor devices.

Theoretical Framework or Literature Review

We consider the change in the properties of the discrete filter 𝐷0 ∙ (𝑗 ∙ 𝜔, 𝑇0) in the frequency domain when the quantization frequency changes N times: fN=f0/N, with constant numerical parameters. It is known that a fairly good approximation of the transcendental FR 𝐷0 ∙ (𝑗 ∙ 𝜔, 𝑇0) = 𝐷0(𝑒𝑗∙𝜔∙𝑇𝑜∗𝑁) in the discrete system in the low frequency range 0 < 𝜔 < 2 ∙ 𝑓0 give pseudo-frequency characteristics (Shamrikov, 1985) constructed as a function of the parameter λ = 2 ∙ 𝑓0 ∙ 𝑡𝑔 𝜔 , where λ is the absolute 2∙𝑓0 pseudo-frequency practically coincides with the real frequency 𝜔, if 𝜔 → 0. We note that the parameter λ in some works on digital signal processing is described by the term "digital frequency". Further we consider the calculation of the parameters of such DF from the condition that the pseudo-frequency characteristics of the reference DF at a given quantization frequency f0 and pseudo-frequency characteristics of such a filter at a sampling frequency f0/N were close within the specified range. Such a criterion can give a good coincidence of the characteristics of a similar filter and a reference filter in the real frequency band 0 <𝜔 < 𝑚𝑖𝑛 (2/𝑇0, 2/𝑇0 ∙ 𝑁). This follows from the fact that the absolute pseudo-frequency λ and real frequency 𝜔 within the specified range are close enough. Using the transition to w-region, where w is the parameter of the bilinear transformation 𝑧 = 1+T0⋅𝑤/2 1−𝑇0⋅w/2 (Shamrikov, 1985), and equating the characteristics of the reference DF at the calculated frequency f0 and a similar DF at a changed frequency fN we will get the following:

whence we write the transfer function of similar a filter in w domain in the following form:

Using the transfer function DN(w) of a similar filter in w domain from (6), and applying the inverse w- transformation for the changed frequency fN, we pass to z-transfer function DN(z) of similar a filter, taking into account the changed frequency f0/N:

Let the reference DF be described by z-transfer function D0(z) of the form (1). Then from expression (6) after transformations we obtain DN(w) and, substituting the formula for the inverse w-transformation 𝑤 = 2⋅(𝑧−1) 𝑇0⋅𝑁 (𝑧+1),we will get the following:

To simplify the expression, we introduce the conversion parameter 𝑅 = 𝑁−1 = 𝑓0−𝑓𝑁, associated with a 𝑁+1 𝑓0+𝑓𝑁 change in the quantization frequency (R=0 at a constant frequency when fN=f0). Taking into account the entered coefficient, dividing the numerator and denominator DN(z) on(𝑁1) 𝑛 ≠ 0 after algebraic transformations, we obtain z - the transfer function of similar a filter in the form of the ratio of polynomials AN(𝑧−1 ) and BN(𝑧−1 ):

where𝑧 = 𝑒𝑗∙𝜔∙𝑇𝑜∗𝑁. Formula (9) defines the way to calculate the coefficients {𝑎𝑁𝑖, 𝑏𝑁𝑗} of z - transfer function of similar DF by coefficients {𝑎, 𝑏} of the reference filter and the recalculation parameter 𝑅 = 𝑁−1 from the condition 𝑖 𝑗 𝑁+1 of closeness of their characteristics at absolute pseudo-frequencies. For an arbitrary order of DF we obtain from expression (9) rather cumbersome relations for determining the coefficients, which at a high order n of the filter will require complex calculations in real time for tuning. In this connection, let us further consider a computational algorithm for obtaining parameters of such a filter: numerator and denominator coefficients (9). From expression (9) we obtain the algorithm for calculating polynomial coefficients AN(𝑧−1) of numerator of the transfer function DN(z) of such a filter, using the representation of the polynomial in the form of Horner's method (Aho et al., 2000):

The polynomial of the denominator BN(z) is similarly obtained for this purpose in formula (10) it is necessary to replace the coefficients 𝑎𝑖 of the numerator polynomial by the coefficients 𝑏𝑖. We denote the polynomials of degree j of 𝑧−1 in square brackets of expression (10) as 𝑀1(𝑧−1, 𝑗), and the polynomials in curly brackets as 𝑀2(𝑧−1 , 𝑗). At the same time, the degree j of the polynomials M1 and M2 varies from 0 to n, and the elements of the polynomials are of degree 𝑧−1 with numerical coefficients. Initial values of polynomials at j=0 are defined as: Subsequent polynomials of degree 1 ton from Horner's method are calculated using the recurrence formula:

Calculation of polynomials 𝑀1(𝑧−1, 𝑗) and 𝑀2(𝑧−1 , 𝑗) ends at the step j=n, and the resulting polynomial 𝑀1(𝑧−1 , 𝑛) corresponds exactly to the numerator 𝐴𝑁(𝑧−1 ) of z-transfer function DN(z) of similar a filter. Similarly we can calculate the polynomial 𝐵𝑁(𝑧−1 )-denominator DN(z), by changing in the algorithm (13) the coefficients of 𝑎𝑖 to 𝑏𝑖⋅ Thus, according to the known polynomials 𝐴(𝑧−1 ) and 𝐵(𝑧−1 ) of z-transfer function D0(z) of the reference filter and the values of the actual quantization period T0*N, or the polynomials can be calculated using the scaling parameter R 𝐴𝑁(𝑧−1 ) and 𝐵𝑁(𝑧−1 ), defining the numerator and denominator of z- transfer function DN(z) of such a filter. Due to the linear relation between the coefficients {𝑎𝑖,𝑏𝑖} of the reference filter and the coefficients {𝑎𝑁𝑖,𝑏𝑁𝑖} of similar a filter, it is convenient to represent this relation in the form of matrix relations:

where A0, B0 are the (n+1)-dimensional parameter vectors {𝑎𝑗}, {𝑏𝑗} of DF reference, AN and BN are the (n+1)-dimensional parameter vectors {𝑎𝑁𝑗}, {𝑏𝑁𝑗 } of a similar DF, Q(R) – (𝑛 + 1) × (𝑛 + 1) is the transformation matrix. It is characteristic that the elements of the matrix Q(R) transformations have, for a given order n, a constant structure and depend only on the parameter R. These elements are polynomials of degree R, they are not higher than n, and the coefficients of these polynomials, integers, are binomial coefficients and depend only on the order n of the filter. Thus, for the second-order filter example (n=2), calculating by formula (10), we obtain:

whence the scaling matrix Q(R) is defined by the following expression:

Methodology

As it was shown above, to ensure the constancy of DF characteristics when the quantization frequency changes, it is necessary in the process of filter operation to carry out the tuning, based on the recalculation of DF parameters. We consider the possible logic of such a tunable filter and the method of its calculation. At the same time, we take into account that the structure of the matrix Q(R) depends only on the order of the digital filter, and the numerical parameters, the elements of the matrix, depend on the specific values of R. Due to the fact that the matrix Q(R) depends only on the degrees of the parameter R, it is convenient to represent this transformation matrix in the form of decomposition:

where 𝑞𝐽 (𝐽 = 0.. 𝑛) – (𝑛+ 1) × (𝑛 + 1) are the matrices, the elements of which are integers, defined by products and sums of binomial coefficients. Taking into account the decomposition (18), the parameters {𝑎𝑁𝑗}, {𝑏𝑁𝑗} of a similar filter from formula (13) is determined by the following expression:

In the formulas (19) 𝑞𝐽 ∙ 𝐴0 and 𝑞𝐽 ∙ 𝐵0 − (𝑛 + 1) – dimensional vectors, their components are known in advance, because they depend on the parameters {𝑎𝑗}, {𝑏𝑗} of the reference DF and elements of constant matrices 𝑞𝐽(𝐽 = 0. . 𝑛). Thus, to recalculate the parameters of DF when changing the quantization period it is enough to 2 ∙ (𝑛 + 1) multiplications of known vectors 𝑞𝐽 ∙ 𝐴0 and 𝑞𝐽 ∙ 𝐵0 (𝐽 = 0. . 𝑛), placed in memory for different degrees of the current recalculation parameter 𝑅𝐽 and 2∙n additions of the obtained products according to the ratios (18). Thus, for the above example of a second-order DF we obtain:

We use the assumption that the period of incoming data can change significantly during operation, but such changes are episodic and the filter basically works in quasistationary mode with a constant step of sampling TN=T0*N and constant parameters. Besides, from practical considerations it is clear that at small enough change of sampling rate | f0-fN | < eps, (eps is a small value) due to the continuity of the parameters does not require. It is assumed, that parameters A0 and B0 of the reference filter are calculated beforehand, and the filter is realized in the form of the program of calculations, similar to (2), for the digital calculator. In this case, the organization of the restructuring is proposed to be carried out at the preliminary stage (labor-intensive operations), and directly in the process (in real time). Pre-calculation stage:

1. to determine the threshold eps value of the frequency change f0-fN, above which the filter parameters must be rebuilt;


2. to calculate the constants (𝑛 + 1) × (𝑛 + 1) - dimensional matrices q0, q1,.., qn, determining the recalculation matrix, and place them in memory, for which: to calculate the recalculation matrix Q(R) from (13) as a function of the parameter R using the computational algorithm from (7), to get matrices qJ - decomposition components Q(R) from (10) and to write them into the memory block – 2 ∙ 𝑛2 ∙ (𝑛 + 1) integers;


3. to write parameters A0, B0 of the reference filter to RAM – 2 ∙ (𝑛 + 1) real numbers;


4. to calculate 2 ∙ (𝑛 + 1), (𝑛 + 1) -dimensional matrices 𝑞𝑗 ∙ 𝐴0 and 𝑞𝑗 ∙ 𝐵0 required for recalculation of parameters by formulas (19,20), place them in the filter memory.

Real-time filter operation:

1. to make episodic measurements of the frequency of incoming information fN;


2. if the frequency change threshold is exceeded: |fN-f0|> еps, to make a quick recalculation of the filter parameters, to do this:

a. to calculate the current recalculation parameters 𝑅 = (𝑓0−𝑓𝑁) , 𝑅2, . . 𝑅𝑛; (𝑓0+𝑓𝑁)

b. to calculate the parameters 𝑎𝑁𝑖, 𝑏𝑁𝑖 of similar filter for specific values of the parameter R and its degrees, using formulas (19,20), and vectors 𝑞𝑗 ∙ 𝐴0 and 𝑞𝑗 ∙ 𝐵0, placed in the filter memory in the preliminary stages;

c. to carry out a realignment and to use the calculated parameters 𝑎𝑁𝑖, 𝑏𝑁𝑖 in the current operation of DF.

Further, the parameters 𝑎𝑁𝑖, 𝑏𝑁𝑖 and the current frequency fN can be used as reference parameters 𝑎𝑖, 𝑏𝑖, and fN is f0 and then we go to the step 3. In the above procedure only items 5 and 6 are performed in real time, i.e. the realignment requires a single execution of uncomplicated and fully defined multiple arithmetic operations of multiplication and addition, while the main time-consuming operations are performed at the calculation stage.

Results

Thus, an algorithm for the operational restructuring of the numerical parameters of the filter based on information about the time intervals of information receipt has been developed. Such a restructuring makes it possible to ensure the constancy of the frequency characteristics of the digital filter with a significant change in the time periods TN=T0*N of information receipt and is implemented by simple computational operations. A technique for implementing the developed recalculation algorithm is proposed, which uses the most cumbersome calculations at the preliminary calculation stage, and during real-time it uses restructuring, performs only single arithmetic calculations with numbers placed in the device's memory. The developed proposals can be used for an arbitrary order of digital filters, while changing the frequency fN of information receipt within a wide range. The approaches proposed in the work for organizing the restructuring of the parameters of DF were tested on a model example of a fourth-order filter. At the same time, preliminary calculations of the conversion matrices and operational calculations in real time were carried out with a limited calculation accuracy corresponding to 8 binary digits. The change in the quantization period T0 corresponded to the range of change in the parameter N=0.3-3, with a corresponding change in the parameter R. Figure 1 shows the frequency characteristics of the fourth-order low-frequency DF that implements the Butterworth distribution (Oppenheim & Schafer, 2018; Jing et al., 2019). Here, the solid thick line defines FR of the reference filter at the calculated frequency f0, the thin line is the response when the quantization frequency fN changes (N=2.5), the dashed line illustrates the response of a similar (tuned) filter. As can be seen from the Figure 1, when changing the step of discreteness of data receipt without restructuring the parameters, the transformation characteristics change significantly. The use of tuning the filter parameters makes it possible to ensure the practical constancy of DF characteristics with a significant change in the quantization frequency, which confirms the validity and prospects of using the developed approaches to ensure the constancy of DF properties when changing the time parameters of the data arrival.

Figure 1. Amplitude-frequency characteristics of DF A possible specific technical implementation of the described approaches is implemented in a digital device according to patent No. 2631976, described in (Belonogov, 2017; Dimc et al., 2017).

Conclusion

Based on the research presented in the work, the following provisions can be formulated. • changing the frequency fN of DF information arrival relative to the calculated frequency f0 with its fixed parameters can significantly worsen the required signal processing properties, and when using the results of DF in a closed loop, significantly degrade the system quality indicators. • technique and a computational procedure have been developed that form the law of rearrangement of the parameters of DF when changing the intervals of information arrival, relative to the calculated values of T0. This procedure is used at the stage of calculating the system, and simplified arithmetic algorithms adjust the parameters in real time only when the frequency fN-f0 changes, exceeding the pre-calculated threshold. • use of tuning the values of DF parameters to ensure the closeness of the filter characteristics to the characteristics of the reference filter operating at the initial frequency f0, allows maintaining the comparative constancy of the filter characteristics with significant changes in the period T0 of data arrival. At the same time, the algorithms for DF and tuning DF parameters are compact and can be easily implemented on the simplest microprocessor devices.