Digital Integrators with an Examples. C Code and Octave Script
1. Abbreviation
2. Introduction
It follows from (2) that the transfer function of the analog ideal integrator:
or in the frequency space of (7), (10):
Then the magnitude and phase responses of the ideal integrator has the form:
Figure 21: Bode Diagram of the Analog Ideal Integrator
Notes:

The integrator is an unstable element. Instability occurs at zero frequency as a result of the integrator property. For example, if a constant signal is applied to the integrator, then output signal of the integrator goes to infinity:
See the Figure 22:
Figure 22: Integrator Response on the Constant Signal 
Next we consider digital approximations of integrators using zTransform (4)
3. Backward Euler Integrator
If the analog prototype (10) is converted using the Backward Euler Method (5), then Backward Euler Integrator is got:
Note. The same result is obtained if the Impulse Invariant Method is used
Moving from the zTransform to the signal samples:
In this case the area is calculated using the rectangle rule. See the Figure 31
Figure 31: Support of the Backward Euler Integrator
Calculate the area in step 10 using the formula (14), substituting all previous steps into it:
Obviously the area of the rectangles is calculated (see the Figure 31):
Magnitude and Phase Response of the Euler Integrator are calculated from (13) using substitution (8):
Figure 32: Magnitude and Phase Resposes of the Backward Euler Integrator
Comparing Figure 21 and Figure 32, it can be seen that the Phase Response of the Backward Euler Integrator is incorrect. This indicates a poor approximation of this digital integrator.
Now let’s calculate the integral using the Backward Euler Integrator. To do this, 5 samples of the sine function will be taken as a test signal:
Figure 33: Input Test Signal: Five Samples of the Sine Function: x(0), x(1), x(2), x(3), x(4)
The sine function:
Let’s choose, for example, the period and the initial phase of the sine:
Then the area of the sine is from 0 to 5Ts:
Let:
fs = 1000Hz => Ts = 0.001s (21)
T = 6Ts = 0.006s (22)
From (20), (21) and (22) the exact area is got :
Now calculate the area using the formula (14). See the DigitalIntegrators.m:
Euler’s Integrator Quality
Half Period Samples Number = 3
exactArea = 1.354866575281059e03
calculationArea = 8.191520442889925e04
relativeMistakeInPercent = 39.5%
The relative error is almost 40%, the mistake is large!
4. Bilinear Integrator
If the analog prototype (10) is converted using the Bilinear Method (6), the Bilinear Integrator is got:
Moving from zTransform to signal samples:
Here the area is calculated using the trapezoid rule. See the Figure 41
Figure 41: Support of the Bilinear Integrator
Calculate the area in step 10 using the formula (25), substituting all the previous steps into it:
Obviously the area of the trapezoids is calculated in this case (see the Figure 41):
The Magnitude and Phase Responses of the Bilinear Integrator is calculated from (24) using substitution (8):
Figure 42: Magnitude and Phase Resposes of the Bilinear Integrator
Note: Magnitude Response has a significant deviation at high frequencies from the ideal integrator, compare with the Figure 21
Now the integral is calculated with the Bilinear Integrator using the test signal (17) – (23) and Figure 33. See also the Digital Integrators.m:
Bilinear Integrator Quality
Half Period Samples Number = 3
exactArea = 1.354866575281059e03
calculationArea = 1.272305937807317e03
relativeMistakeInPercent = 6.1%
The relative error is approximately 6%.
5. Simpson’s Integrator
To approximate the calculation of an integral, the area under the graph shown in Figure 51, Simpson’s formula (28) can be used:
Figure 51: Support of the Simpson’s Formula
The proof of the Simpson’s formula was given in the article: “The RungeKutta and Hamming Methods for Numerical Solution of Ordinary Differential Equations with an Example. C Code and Octave Script”
Now using the Simpson’s formula, the area is calculated under the graph x(t) in the range t0…t10 shown in Figure 52
Figure 52: Support of the Simpson’s Integrator
The area calculation formula can also be written in recursive form, that is in the form of Simpson’s IIR filter:
zTransform of the Simpson’s Filter:
The Magnitude and Phase Responses of the Simpson’s Integrator is calculated from (31) using substitution (8):
Figure 53: Magnitude and Phase Resposes of the Simpson’s Integrator
The Magnitude Response has an unexpected rise in the upper frequencies, which does not correspond to the ideal integrator, see the Figure 21. This is due to the fact that the integrator is unstable at the fs/2 frequency. It can be seen from (31) that the transfer function has two poles on the unit circle:
which corresponds to the frequencies f1 = 0 and f2= fs/2. Instability at zero frequency is the norm for an integrator and was discussed above. Instability at the fs/2 frequency is a problem and can lead to incorrect integral calculations. Therefore it is important that there are no high frequencies in the input signal spectrum. If a short signal with a frequency of fs/2 is applied to this integrator input then continuous oscillations will occur:
(a) Short input signal with a frequency of fs/2
(b) Continuous oscillations at the output of the Simpson’s integrator
Figure 54: Simpson’s Integrator: Instability at the fs/2 frequency
Now let’s calculate the integral with Simpson’s Integrator, using the test signal (17) – (23) in the Figure 33. See also the Digital Integrators.m:
Simpson’s Integrator Quality
Half Period Samples Number = 3
exactArea = 1.354866575281059e03
calculationArea = 1.394305321397539e03
relativeMistakeInPercent = 2.9%
The relative error is approximately 3%.
6. Tick’s Integrator
Leo Thick’s formula is given and proved in [1], which approximately calculates an integral:
zTransform of the Tick’s Integrator:
The Magnitude and Phase Responses of the Tick’s Integrator is calculated from (34) using substitution (8):
Figure 61: Magnitude and Phase Resposes of the Tick’s Integrator
Note. Just like Simpson’s Integrator, the Tick’s Integrator is unstable at the fs/2 frequency. See the Figure 54.
Now calculate the integral with the Tick’s Integrator using the test signal (17) – (23) in the Figure 33. See also the DigitalIntegrators.m:
Tick’s Integrator Quality
Half Period Samples Number = 3
exactArea = 1.354866575281059e03
calculationArea = 1.375956614105570e03
relativeMistakeInPercent = 1.6%
The relative error is approximately 1.6%. That is the best result among the listed integrators!
7. Conclusions

Integrators Quality can be analyzed and compared using the method from [1]. To do this, the relation will be considered:
Then using (35) and (36):
Figure 71 shows all four graphs together from (37) – (40):
Figure 71: Integrators Quality Compare
Notes:
– The Backward Euler Integrator has not only amplitude deviations from the ideal integrator, but also a phase approximation error (see the Figure 32). The Figure 71 shows only amplitude deviations from the ideal integrator. The Backward Euler Integrator is the worst approximation. Despite this the Backward Euler Integrator is used only in simple Digital Signal Processing (DSP) applications, for example, in PI Regulators
– From the Figure 71 the best integrator is Tick’s Integrator, followed by Simpson’s Integrator, and Bilinear Integrator is in third place

Simpson’s and Tick’s Integrators have additional instability at the fs/2 frequency, which makes them difficult to use in DSP. These integrators can be used, for example, in conjunction with a Low Pass Filter, which removes the upper unstable part of the spectrum.

The integrator has instability at zero frequency (see the Figure 22). If a const signal is applied to the integrator input, then the signal goes to infinity at the output of the ideal integrator. In some applications this effect be had to deal and to use the AntiWindup Method. It consists in the fact that the magnitude of the output signal and previous output signal value(s) are limited to the operating range and does not allow the system to enter saturation
Example of the AntiWindup for Bilinear Integrator:

The quality of integration also depends on the input signal. For example, for a stepwise signal, the best integration result can be got using the Back Euler Integrator, despite the fact that this is the worst approximation of the integrator

I use Bilinear Integrators in the article “Control of Non Linear Object Using the Vostrikov’s Localization Method. C Code and Octave Script”

I advise you to read more on the topic of digital integrators in [1]
8. Octave GNU file DigitalIntegrators.m
The .m script generates the digital integrator solutions and plots. To do this, the following functions are defined in this file:
% Plotting the Magnitude Response of the Second Order Analog Filter: (b1*s^2+b2*s+b3)/(a1*s^2+a2*s+a3)
% Input parameters:
% SB – numerator of the second order section
% SA – denominator of the second order section
% FhighInHz – Frequency space in Hz
function AnalogMagnitudeResponse(SB, SA, FhighInHz, Text)
% Plotting the Phase Response of the Second Order Analog Filter: (b1*s^2+b2*s+b3)/(a1*s^2+a2*s+a3)
% Input parameters:
% SB – numerator of the second order section
% SA – denominator of the second order section
% FhighInHz – Frequency space in Hz
function AnalogPhaseResponse(SB, SA, FhighInHz, Text)
% Plotting the Magnitude Response of the Digital Second Order filter: (b1+b2*z^1+b3*z^2)/(a1+a2*z^1+a3*z^2)
% Input parameters:
% ZB – numerator of the second order section
% ZA – denominator of the second order section
% Tsample – sampling time in seconds
function dig_func_dmr = plot_DigitalMagnitudeResponse(ZB, ZA, Tsample, Text)
% Plotting the Phase Response of the Digital Second Order Filter: (b1+b2*z^1+b3*z^2)/(a1+a2*z^1+a3*z^2)
% Input parameters:
% ZB – numerator of the second order section
% ZA – denominator of the second order section
% Tsample – sampling time in seconds
function dig_func_angle = plot_DigitalPhaseResponse(ZB, ZA, Tsample, Text)
% Plotting the Ratio: abs(DigitalIntegrator/IdealIntegrator) = w*abs(DigitalIntegrator) = 2*pi*f*abs([(b1+b2*z^1+b3*z^2)/(a1+a2*z^1+a3*z^2)])
% Input parameters:
% ZB – numerator of the second order digital integrator
% ZA – denominator of the second order digital integrator
% Tsample – sampling time in seconds
function dig_func_ratio = plot_DigitalIntgratorQuality(ZB, ZA, Tsample)
% Float Table to File
% Support CCode
% Input parameters:
% ParameterVector – data array
% FileNameString – output file name
function FloatParamVector2file(ParameterVector, FileNameString)
9. C Language DigitalIntegrators.c/h
The demo SW is realized the digital integrator solutions. See the files DigitalIntegrators.c/h
Functions:
/*
IIR Filter with Transposed Form II
Input
float input – filter input
Return
float – filter output
*/
float iir_filter(float input, float const* all_coef_ptr, float* all_history_ptr, uint8 sos_nmb)
/*
Initialization of the Integrators
Clear the working arrays
Input: void
Return: void
*/
void IntegratorsInit(void)
/*
Euler’s Integrator
Input
float input – input value
Return
float – output value
*/
float EulerIntegrator(float input)
/*
Bilinear Integrator
Input
float input – input value
Return
float – output value
*/
float BilinearIntegrator(float input)
/*
Simpson’s Integrator
Input
float input – input value
Return
float – output value
*/
float SimpsonIntegrator(float input)
/*
Tick’s Integrator
Input
float input – input value
Return
float – output value
*/
float TickIntegrator(float input)
10. Download the DigitalIntegrators.c/h/m
You can download the files:
DigitalIntegrators.m
DigitalIntegrators.c
DigitalIntegrators.h
with the button: