Digital Integrators with an Examples. C Code and Octave Script

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 2-1: 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 2-2:

  
  Figure 2-2: Integrator Response on the Constant Signal

• Next we consider digital approximations of integrators using z-Transform (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 z-Transform to the signal samples:

In this case the area is calculated using the rectangle rule. See the Figure 3-1

Figure 3-1: 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 3-1):

Magnitude and Phase Response of the Euler Integrator are calculated from (13) using substitution (8):

Figure 3-2: Magnitude and Phase Resposes of the Backward Euler Integrator

Comparing Figure 2-1 and Figure 3-2, 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 3-3: 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.354866575281059e-03
calculationArea = 8.191520442889925e-04
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 z-Transform to signal samples:

Here the area is calculated using the trapezoid rule. See the Figure 4-1

Figure 4-1: 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 4-1):

The Magnitude and Phase Responses of the Bilinear Integrator is calculated from (24) using substitution (8):

Figure 4-2: 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 2-1

Now the integral is calculated with the Bilinear Integrator using the test signal (17) – (23) and Figure 3-3. See also the Digital Integrators.m:

Bilinear Integrator Quality
Half Period Samples Number = 3
exactArea = 1.354866575281059e-03
calculationArea = 1.272305937807317e-03
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 5-1, Simpson’s formula (28) can be used:

Figure 5-1: Support of the Simpson’s Formula

The proof of the Simpson’s formula was given in the article: “The Runge-Kutta 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 5-2

Figure 5-2: 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:

z-Transform of the Simpson’s Filter:

The Magnitude and Phase Responses of the Simpson’s Integrator is calculated from (31) using substitution (8):

Figure 5-3: 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 2-1. 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 5-4: 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 3-3. See also the Digital Integrators.m:

Simpson’s Integrator Quality
Half Period Samples Number = 3
exactArea = 1.354866575281059e-03
calculationArea = 1.394305321397539e-03
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:

z-Transform of the Tick’s Integrator:

The Magnitude and Phase Responses of the Tick’s Integrator is calculated from (34) using substitution (8):

Figure 6-1: 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 5-4.

Now calculate the integral with the Tick’s Integrator using the test signal (17) – (23) in the Figure 3-3. See also the DigitalIntegrators.m:

Tick’s Integrator Quality
Half Period Samples Number = 3
exactArea = 1.354866575281059e-03
calculationArea = 1.375956614105570e-03
relativeMistakeInPercent = 1.6%

The relative error is approximately 1.6%. That is the best result among the listed integrators!

 

 

7. Conclusions

1. 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 7-1 shows all four graphs together from (37) – (40):

   

   Figure 7-1: 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 3-2). The Figure 7-1 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 7-1 the best integrator is Tick’s Integrator, followed by Simpson’s Integrator, and Bilinear Integrator is in third place
   

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

3. The integrator has instability at zero frequency (see the Figure 2-2). 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 Anti-Windup 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 Anti-Windup for Bilinear Integrator:

   
   

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

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

6. 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 C-Code
% 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:

 

 

11. Literature / References

[1] R. W. Hamming “Digital Filters”, Third Edition, Dover Publications, INC. Mineola, New York, 1997, ISBN 10: 0-486-65088-X / ISBN 13: 978-0486650883