Design of Second Order Recursive Digital Integrators with Matching Phase and Magnitude Response

Location of poles and zeroes greatly affect phase response and magnitude response of a system. Recently, pole-zero optimization emerged as an effective approach to approximately match magnitude response of a system with that of an ideal one. In this brief, a methodology for the design of linear phase integrators and ones with constant phase of −90 degree is proposed. The aim of this method is to simultaneously attain multiple objectives of magnitude and phase optimization. In this method, magnitude response error is minimized under the constraint that the maximum passband phase-response error is below a prescribed level. Examples are included to illustrate the proposed design technique.


Introduction
An integrator is a system whose output signal is the time integral of its input signal.It can be mathematically modelled as, H i (ω) = 1/jω where j = √ −1 and ω is the angular frequency in radians per second.Integrators find immense applications in the areas of signal processing, bio-medical engineering, radar engineering, sonar engineering, control system etc.Multitude of techniques have been developed for the design of recursive and non-recursive digital integrators.It is interesting to note that for linear phase integrators, Maximum Percentage Relative Error (MPRE) in magnitude response has rolled down from 5 % to 0.48 % in the decade 2005 − 2015.
Several designs of digital integrator have been proposed in the literature [1][2][3][4][5][6][7].Among these, the simplest approximations to the desired frequency response are Rectangular, Trapezoidal and Simpson digital integrators.It can be easily noticed that, as ω → 0, H (e jω ) → 1/jω for these approximations.These digital integrators can be sufficient for the integration of oversampled signals (signals sampled much above the Nyquist rate), but a deeper investigation suggests that there is still a possibility of better designs for signals that are critically sampled at the Nyquist rate.
Digital integrators have been designed using quadrature rules, such as the Newton-Cotes and Gauss-Legendre rules by Ngo and Tseng [3][4][5][6].These methods, however, are complex to design and implement due to the usage of fractional sampling rates, hence lacking computational efficiency.Lagrange interpolators have been suggested to elevate some of these problems.Most recently Tseng et al. have proposed to implement the fractional delays in the Hartley transform domain [7].
One of the design methods which has become very popular among researchers is Iterative Optimization Method [8][9][10][11][12][13][14][15][16][17][18][19][20].It is widely used to improve the performance of a system by reducing its runtime, bandwidth, memory requirement, or some other property.Optimization methods such as Linear Programming, Simulated Annealing, Genetic Algorithm, and Pole-Zero optimization have been used earlier to design Infinite Impulse Response (IIR) digital integrators and differentiators.Papamarkos-Chamzas [8] have used Linear Programming optimization method to design digital integrators.Al-Alaoui [9] has also proposed a family of digital integrators by using interpolation and Simulated Annealing optimization method.Upadhyay-Singh (US) have proposed recursive wideband digital integrator for 0.48 % MPREs in magnitude responses over almost the full Nyquist band except near to ω = π [10].Genetic Algorithm has been exploited in [16] to obtain a class of second order linear phase integrators.In [17], Jain-Gupta-Jain have used Minimax and Pole, Zero, and Constant Optimization Methods to obtain second, third and fourth order IIR digital integrators.Later, in [18], Gupta et al. (GET) have proposed recursive wideband digital integrators using Modified Particle Swarm Optimization (MPSO).In [19], Jalloul and Al-Alaoui have employed Particle Swarm Optimization to propose designs of integrators and differentiators.Upadhyay, in [20], has used Pole-Zero Optimization to achieve integrators with relative error in magnitude response not exceeding 0.37 %.In [21], efficient design of FIR digital differentiator using the L 1optimality criterion is proposed.
Often application in controls, wave shaping, oscillators and communication require a constant 90°phase for differentiators and −90°phase for integrators.In [22], Al-Alaoui cascaded differentiator and integrator operator with fractional advance and delay respectively to obtain constant 90°phase for differentiators and constant −90°phase for integrators.Al-Alaoui also showed that doubling the sampling rate improves the magnitude response.Combining the two actions improves both the magnitude and phase responses.It should be noted that the approach applies to other differentiators and integrators with linear phases, or approximately linear phases, such as the second-order Al-Alaoui integrators and differentiators.However, the approach is limited only to linear phase differentiator and integrator.Also, fractional delays and advances are complex to implement.In this brief, our goal is to improve the performance and computational efficiency of digital integrators as a standalone system.
Iterative constrained optimization has already been proven a superior method of obtaining close-to-perfect magnitude responses.However, cost function defined in these methods only considers magnitude error with no account of phase error.For real-time applications, further reduction in magnitude error is not as important as reduction in phase error is.
Every recursive digital filter can be completely specified by its poles and zeros with suitable scaling.Poles and zeros give useful insights into a filter's frequency response, and can be used as the basis for digital filter design.This paper attempts at developing a design method which governs constraints on the location of poles and zeroes in order to achieve multi-objectives of phase and magnitude responses, both together.

Problem Formulation
The paper addresses two design objectives.Firstly, the design of stable recursive linear phase digital integrator and secondly, the design of stable recursive digital integrator with a constant phase of −90°.These designs, however, involve non-converging objectives of magnitude and phase optimization and there exist a trade-off among the different objectives.The trade-off parameters considered in the designs are: passband magnitude error and passband phase deviation, passband cut-off frequency in phase response (ω p ): which is defined as range of frequency for which phase deviation ≤ δ.
Design objectives for linear phase integrators can hence be summarized as: 1. Low wideband magnitude error 2. Low passband phase deviation (constant group delay) Design objectives for integrators with constant phase of −90°a re given below: 1. Low wideband magnitude error 2. High Passband cut-off frequency in phase response

Low passband phase deviation
Let, H l (z) denotes the transfer function of linear phase integrator and H c (z) denotes the transfer function of integrator with constant phase of −90°.
A generalized digital recursive transfer function of second order is of the form, given in (1).
where z 1 , z 2 are zeroes and p 1 , p 2 are poles of the considered transfer function.Here, I o is a multiplier constant or scaling factor.The frequency response I (e jω ) of the considered system can be written as in (2), where Error function to be minimized is given in (3).
The main reason for expressing the magnitude constraints in terms of magnitude-square instead of magnitude is to avoid the appearance of a denominator for the gradients of the magnitude constraints, which the use of the simple magnitude would have made unavoidable.
A necessary and sufficient condition for the stability of a causal Linear Time Invariant (LTI) digital IIR filter is that all poles of its irreducible transfer function lie strictly inside the unit circle.

Solution Methodology
Consider a polynomial equation x = (z − a), where a ∈ R. Substituting z = e jω , we get phase of x, given in (4).location of root Groupdelay of x is defined as: From (5a), it is noted that: Therefore, indicating that greater the value of a, better the phase linearity of the system.
We next require to compute values of φ x (ω) at extreme frequency points in the Nyquist Band (i.e. at ω = 0 and at ω = π) for different locations of root.Table 1 summarizes values of φ x (ω) for different root locations.

Design of Linear Phase Integrators
Figure 1 illustrates the phase effects of zeroes and poles to obtain linear phase characteristics for a discrete time system.It is observed from this figure that in order to have a linear phase response for a system with transfer function H l (z), phase function of one of the poles p 1 , i.e, φ 3 (ω) should be equal to 90°at ω = 0 and decrease monotonically henceforth with constant slope for the complete range of Nyquist frequency.Another pole p 2 should be located strictly inside the unit circle to satisfy stability constraints for the system.One of the zeroes z 2 should be so chosen that it neutralizes the phase effects of p 2 .The phase response of zero z 1 , i.e, φ 1 (ω) is required to be monotonically decreasing, linear function of frequency with slope much greater than that of φ 3 (ω), so that, when phase characteristic of p 1 is subtracted from that of z 1 , i.e, φ 1 (ω) − φ 3 (ω), linear phase characteristic is obtained for the overall transfer function H l (z).
Based on this, required phase response of p 1 can be mathematically modelled by the following conditions: From Tab. 1, it is observed that for pole location z = 1, first two of the above mentioned conditions are fulfilled.
Hence, it is concluded that pole p 1 is located at z = 1.Isolated poles on the unit circle may be called marginally stable.
The impulse response corresponding to a single pole on the unit circle never decays, but neither does it grow.
Next, constraints on phase function of z 1 are described below: For the first condition to hold true, zero z 1 can take any value in the interval (−∞, −1), as concluded from Tab. 1. Assume z 1 lies in the interval (−α, −1).Solution for α can be obtained by solving the equation: κ = 1 results into a perfect linear phase response while κ = 0 leads to a constant phase response of −90°.
Different solutions of α can be obtained for different defined values of κ.
An LTI filter is stable if and only if its poles are strictly inside the unit circle (|z| = 1 ) in the complex z plane.In particular, a pole p outside the unit circle (|p| > 1 ) gives rise to an impulse-response component proportional to p n which grows exponentially over time n.Therefore, pole p 2 must lie in the interval (-1,1).To achieve linear phase characteristics for the system, p 2 and z 2 should lie close enough to each other.Let ∆p be the distance between the two.For faithful neutralization of the two phase effects, |∆p| should not exceed 0.1 In conclusion, to obtain linear phase characteristics for a second order system, p 1 should be located at z = 1.p 2 should lie in the interval (-1,1).z 2 = p 2 + |∆p|, where |∆p| does not exceed 0.1.z 1 should lie in the interval (−α, −1) where α is linearly related to user-defined parameter κ.

Design of Integrators with Constant Phase of −90°N
ow that we have estimated location of poles and zeroes of H l (z) to obtain linear phase characteristics for the system, we would now like to investigate location of poles and zeroes to achieve a constant phase response of −90°for a system with transfer function H c (z).Let z 1 , z 2 , p 1 and p 2 denote the zeroes and poles of H c (z).
From Fig. 1 It is to be noted that phase function of pole p 1 suffers the same constraints as p 1 and it continues to be located at z = 1.Hence, ∆ 3 = 0. Now, summation of phase shifts caused by z 1 , z 2 and p 2 should be equal to 26.56°, given by (7), where ∆ is tolerable angular phase shift.At first, consider |∆ 2 | and ∆ small enough to be neglected.So, (7) becomes, In (8a), z 1 ≥ −1.So, maximum value of ∆ 1 is achieved when Relationship between root location a and corresponding maximum angular phase deviation (∆) caused by it is derived as follows: Maximum Angular Phase Deviation (∆) can be defined as slope of a tangent to the phase response curve at ω = 0.By above definition, Similarly, relationship between root location a and corresponding maximum deviation caused (δ) by it is derived as follows: Maximum Phase Deviation (δ) can be defined as the value of phase obtained at a stationary point on the phase curve at which the tangent changes from a positive value on the left of this point to a negative value on the right.
Solving for ω in the equation dφ x (ω) dω = 0 we get

By definition
This gives us Solution of ( 9) and (10b) gives relationship between ∆ and δ.
For different values of δ (and hence ∆) chosen, different solution sets for ( β, γ) can be derived.

Design Approach
The typical approach in designing of a digital IIR integrator is to minimize the maximum amplitude-response error and maximum phase-response error.The optimization can be carried out by minimizing the magnitude error under the constraint that the passband phase error δ is within prescribed level.
Iterative optimization methods are often the only choice for non linear equations.It is a mathematical procedure that generates a sequence of improving approximate solutions until convergence is reached.Genetic Algorithm (GA) is already used in the literature for designing digital differentiators and integrators [13].In this work, GA is used in a similar fashion as in [13].However, the optimization works on optimizing the gain factor and the locations of poles and zeroes of the digital filter as opposed to optimizing the coefficients of the numerator and the denominator.
Steps followed in the design of H l (z) and H c (z) are given below: Step 1: For the design of H l (z), define a value of κ.Then, solution of ( 6) gives the value of α.
Step 2: Run the iterative optimization with poles and zeroes set to vary within the respective prescribed intervals for both the designs, as discussed in previous section.Cost function is given in (3).The optimization runs until convergence is reached.
The constraints on the location of poles and zeroes ensure that the phase-response error is under a prescribed level.The frequency range is defined as 0 ≤ ω ≤ π radians/second.

Design Examples
In real world applications, most of the optimization problems involve more than one objective to be optimized.The objectives in most of engineering problems are often conflicting.In the case, one extreme solution would not satisfy both objective functions and the optimal solution of one objective will not necessarily be the best solution for other objective(s).Therefore, different solutions will produce trade-off between different objectives and a set of solutions is required to represent the optimal solutions of all objectives.Several design examples are included which demonstrate the effectiveness of the design technique.

Linear Phase Integrators
In this section, three designs of linear phase digital integrators have been considered for different values of κ.Figures 2 and 3 illustrate the magnitude and phase responses of proposed integrators H l1−l3 (z).
The design out performs other linear phase design examples proposed in this section in terms of phase linearity over the full Nyquist band and has reasonably well magnitude response for complete Nyquist frequency range.
Example 3: κ = 0.9 Solving (6) for κ = 0.9, α comes out to be 10.Proposed integrator is given in (15).The design outperforms other design examples in terms of magnitude response and also has considerably good phase response over wideband.

Integrators with Constant Phase of −90°T
his section demonstrates three designs of integrators with constant phase of −90°for different values of tolerable maximum phase deviation δ.Figures 4 and 5 show the magnitude and phase responses of designed integrator H c1−c3 (z).
Phase deviation for this design does not exceed 2°in the frequency range 0 ≤ ω ≤ 0.88π.So, cut-off frequency in phase response (ω p ) turns out to be 0.88π.The design, however, can not be considered magnitude efficient with Maximum Magnitude Relative Error (MMRE) of 0.5025 over wideband.
The design has considerably accurate magnitude response within δ = 5°in phase response in the frequency range 0 ≤ ω ≤ 0.69π.Here, ω p comes out to be 0.69π, making it attractive for real-time applications.
The design excels in magnitude response and has δ = 10°in phase response in the frequency range 0 ≤ ω ≤ 0.76π.ω p for this design is 0.76π Though, it is disappointing that none of the design example achieves a constant phase response of −90°over full Nyquist band, it is worth noticing that passband cut-off frequency ω p or bandwidth of the system largely depends on location of zero z 1 .The dependence is depicted below: For a given location of z 1 , approximate passband for the system could be predicted.Conversely, location of z 1 could be determined for a prescribed value of ω p .

Performance Measure
SA is another popular optimization technique that is widely used in the literature for obtaining designs for integrators and differentiators [9].In this work, SA is implemented to simulate integrators proposed in this section.H l1 (z) and H c2−c3 (z) are used as initial guess in the optimization algorithm.
Designs of digital integrators proposed using SA are of the general form, as given in (19).
Proposed linear phase integrators are tabulated in Tab. 2 and ones with constant phase of −90°are tabulated in Tab. 3.
Figures 6 and 7 show the magnitude and phase responses of proposed linear phase integrators H l3−l7 (z), H JGJ (z), H 2MO+PZC (z), H get (z), H up (z) digital integrator with the   ideal one for the complete Nyquist frequency range, respectively.Further, Figures 8 and 9 show the magnitude and phase responses of proposed integrators with constant phase of −90°H c2−c6 (z), Ngo [3] and Tseng [5] digital integrator with the ideal one for the complete Nyquist frequency range, respectively.For proposed and existing linear phase integrators, δ and MMRE, over full Nyquist band, are given in Tab. 4 and that for proposed and existing integrators with constant phase of −90°, for 0 ≤ ω ≤ 0.71π and 0 ≤ ω ≤ 0.76π are given in Tab. 5.

Integrators
Max δ(ω c = 0.71π) Max δ(ω c = 0.76π) Max MRE H c2 (z) H l3 (z) has an MMRE of 0.0273 over wideband, better than the existing or proposed integrator designs.It exhibits a maximum phase deviation of just 3.6089°over complete Nyquist range.H l4 (z) exhibits an MMRE of 0.0433 and a maximum phase deviation of 3.4775°.H l5 (z) and H l6 (z) have an MMRE < 0.02 for 0 ≤ ω ≤ 0.8π and an MMRE of 0.0341 and 0.0398 over wideband, respectively.They exhibit a maximum phase deviation of 3.9078°and 3.7995°o ver complete Nyquist band, respectively.With an MMRE of 0.0580 over wideband, H l7 (z) has relatively low magnitude error for 0.7π ≤ ω ≤ 0.95π.It exhibits maximum phase deviation of 3.8579°over wideband and outperforms all proposed and existing integrator designs for 0.7π ≤ ω ≤ π in phase response.
With reference to Figures 12 and 13, H Ngo (z) and H Tseng (z) perfectly approximate ideal phase response of −90°f or ω upto 0.3π after which their phase deviation drastically increases for the rest of the Nyquist band.H Tseng (z) exhibits excellent magnitude response upto ω = 0.3π but then fails to maintain this for the rest of the Nyquist frequency range.H Ngo (z), however, has an MMRE of 0.0643 over wideband.All proposed designs H c2−c6 (z) have MMRE better than that of H Tseng (z) over wideband.
H c4 (z) outperforms existing and proposed integrators in magnitude response with an MMRE of just 0.0279 over wideband.It is followed by H c3 (z) which has a reason-ably well magnitude response for full Nyquist band with an MMRE of 0.0367.These designs, however, do not stand best in phase response but manage to maintain δ < 10°for 0 ≤ ω ≤ 0.77π and 0 ≤ ω ≤ 0.76π respectively.These designs could be termed as magnitude-dominant designs.
H c5 (z) and H c2 (z) out stand existing and other proposed integrator designs in phase response for 0 ≤ ω ≤ 0.71π with δ < 5°up to ω = 0.71π and δ approximately equal to 5°upto ω = 0.70π, respectively.H c5 (z) outperforms H Ngo (z) for 0 ≤ ω ≤ 0.69π in magnitude response and exhibits an MMRE of 0.1082 over wideband.H c2 (z) lacks in magnitude response with an MMRE of 0.1460 for complete Nyquist band.These designs could be called phase optimal designs.Amongst proposed integrator designs, H c6 (z) exhibits worst magnitude response over wideband and has an MMRE less than 0.1 up to ω = 0.8π and an MMRE of 0.2514 over full Nyquist band.But it is interesting to note that it has δ < 5°for 0 ≤ ω ≤ 0.7π and δ < 10°for 0 ≤ ω ≤ 0.76π.This can be considered as best trade off between δ and ω c among all other proposed designs which either have δ < 5°u pto ω = 0.7π and δ < 10°for ω < 0.76π or have δ < 10°f or ω = 0.76π and δ < 5°for ω < 0.7π.It exhibits a unique bandwidth optimality which none of the other design does.
Extensive simulations have been carried out to verify the system's response of proposed integrators but the actual characteristics of the system, in hardware implementation, is governed by the bit-resolution of the processor.The systems coefficients are required to be converted in binary for processing.The coefficients which are multiple of 5 could be represented in binary by finite number of bits while others require infinite number of bits for their binary representation, hence undergoing digital word length effect, owing to which, systems response might worsen.Now a days, 64-bit processors are commonly used for real-time applications.In such a scenario, impact of digital word length effect of coefficients on system's response will be minimal.The poles of all the proposed integrators have been constraint to be located at z = 1, which ensures the response of the proposed integrators match the ideal one near the origin (ω = 0).

Conclusion
Methodology for the design of second order linear phase integrators and ones with constant phase of −90°is presented.
Proposed integrators H l3−l7 (z) show an excellent phase linearity with δ < 4°over full Nyquist band while simultaneously maintaining low magnitude error.Constant phase of −90°is achieved for proposed integrators H c2−c6 (z) for 0 ≤ ω ≤ 0.73π and they exhibit considerably low magnitude error over wideband.Trade-off among magnitude error, phase deviation and passband cut-off frequency is analysed and results show attractive real-time signal processing applications of proposed integrators as per the requirement of accuracy and bandwidth.
The performances of the proposed integrators are compared with few other state of the art integrators and proven of considerable high accuracy.As a future scope of the research, impact of digital word length of the coefficients on the system's response of proposed integrators could be studied further.

Fig. 1 .
Fig. 1.Illustration of phase effects of zeroes and poles to obtain linear phase characteristics for a system.

Fig. 6 .
Fig. 6.Magnitude responses of proposed and the existing linear phase integrators with the ideal integrator.

Fig. 7 .
Fig. 7. Phase responses of proposed and the existing linear phase integrators with the ideal integrator.

Fig. 8 .
Fig. 8. Magnitude responses of proposed and the existing −90°p hase integrators with the ideal integrator.

Fig. 9 .
Fig. 9. Phase responses of proposed and the existing −90°phase integrators with the ideal integrator.

Fig. 10 .
Fig. 10.Relative Magnitude Errors of proposed and the existing linear phase integrators.

Fig. 11 .
Fig. 11.Relative Phase Errors of proposed and the existing linear phase integrators.
, it is observed that a positive angular shift of 26.56°in phase response of H l (z) causes its linear phase characteristics to transform into a constant phase characteristics of −90°.It can henceforth be argued that cumulative sum of angular phase shifts contributed by poles and zeroes of H c (z) should be equal to 26.56°.
, respectively.Coefficients of proposed integrators with constant phase of -90 degree obtained using SA.
Max Phase Deviation and Max MRE of Integrators with constant phase of -90 degree.Figures 10 and 11show plots of the variation of magnitude error function and phase error function with frequency for the proposed linear phase integratorsH l3−l7 (z), H JGJ (z), H 2MO+PZC (z), H get (z), H up (z) digital integratorsFigures12 and 13show plots of the variation of magnitude error function and phase error function with frequency for integrators with constant phase of −90°designed using our method H c2−c6 (z), Ngo and Tseng integrator.From Figures10 and 11, it is observed that all the proposed linear phase integrator designs H l3−l7 (z) have an MMRE not more than 0.06 and maximum phase deviation less than 4°over wideband.All these designs clearly outperform any of the existing integrators in phase response over complete Nyquist range.H JGJ (z), H 2MO+PZC (z), H get (z) and H up (z) have a maximum phase deviation of 38.4682°, 45.0071°, 37.6768°and 10.4348°for 0 ≤ ω ≤ π and an MMRE of 0.2153, 0.2909, 0.2082 and 0.0149, respectively.