Integrate Routine

Summary

Integrate signal.

Declaration

Procedure Integrate(Src, Dst: TVec; var State: TIntegrateState; Dt: TSample);

Description

Use Simpson's formula to integrate the Src and place the result in Dst. Src can be real or complex. dT defines the sampling period (dT = 1/FS). The State variable holds the initial conditions.


                    j
    Dst[j]  = 1/6* Sum ( Src[i-2] + Src[i-1]*4 + Src[i])*dT)
                   i=0

    j = 0,1...n-1

The integrate routine does not preserve linear phase. An analytical solution for a linear phase integrator does not exists. In general there are two types of applications for integration:

1. The signal has a more or less fixed mean value. (DC offset). An example of such a signal is the signal comming from the accelerometer measuring vibrations. The numerical integration will work well on signals whose mean value is exactly zero (otherwise it will rise or fall in infinity). Because this is often not the case, the signal must be passed through a DC filter first. The DC filter can be IIR or FIR type. An alternative to numerical integration and the Integrate routine is a linear phase integration filter designed with the RemezImpulse routine. Linear phase integrator also removes the DC offset.

2. The signal does not have a mean value. An example of such a signal is the signal comming from the accelerometer measuring the acceleration and deceleration of a driving car. In this case the Integrator routine can be used directly to obtain speed and/or distance from the acceleration data.

See Also
Differentiate
DCFilter
RemezImpulse
KaiserImpulse
FIRImpulse

Example 1

Single block and streaming application example:


uses MtxExpr, Math387, MtxVec, SignalUtils, MtxVecTee, MtxVecEdit;

procedure TForm1.Button1Click(Sender: TObject);
var b,c,h: Vector;
    n,i: integer;
    FS: TSample;
    State2: TIntegrateState;
    State: TDiffState;
begin
    h := Ramp(30,0,1);
    b.Size(h);
    c.Size(h);      
    FS := 1; //sampling frequency

//single block

    FillChar(State,sizeOf(State),0); //initial conditions to zero
    FillChar(State2,sizeOf(State2),0);
    Integrate(h, b, State2, 1/FS);
    Differentiate(b, c, State, 1/FS);
    Drawit(c);

//streaming

    FillChar(State,sizeOf(State),0);  //initial conditions to zero
    FillChar(State2,sizeOf(State2),0);
    n := h.Length div 10;
    for i := 0 to 10-1 do
    begin
        h.SetSubRange(i*n,n);
        b.SetSubRange(i*n,n);
        c.SetSubRange(i*n,n);
        Integrate(h,b, State2, 1/FS);

//            Should be:  b = [0 , 1 , 3, 6, 10, 15, 21,... ]
//            But becomes: b = [0, 0.1666, 1.1666, 3.166, 6.166, 10.166, 15.166, 21.166,... ]
//            because of Simpson

        Differentiate(b,c,State,1/FS);

//            Should be: c = [0,1,2,3,4,5,6....]
//            But becomes: c = [0, 0.08333, 0.5833, 1.5, 2.5, 3.5, 4.5....]

    end;
    c.SetFullRange;
    DrawIt(c);
end;


    #include "MtxVecCPP.h" //MtxVecCPP.cpp must be included in the project
    #include "MtxVecEdit.hpp"
    #include "MtxVecTee.hpp"
    #include "SignalUtils.hpp"
    #include <string.h>

    void __fastcall TForm41::BitBtn1Click(TObject *Sender)
    {
      int n,i;
      double FS;
      TIntegrateState State2;
      TDiffState State;
      Vector b,c,h;

      h = Ramp(30,0,1);
      b->Size(h);
      c->Size(h);
      FS = 1; //sampling frequency

    //single block
      memset(&State,0,sizeof(State)); //initial conditions to zero
      memset(&State2,0,sizeof(State2));

      Integrate(h, b, State2, 1/FS);
      Differentiate(b, c, State, 1/FS);
		
      DrawIt(c);

    //streaming

      memset(&State,0,sizeof(State)); //initial conditions to zero
      memset(&State2,0,sizeof(State2));

    	n = h->Length/10;
      for (i = 0; i < 10; i++)
      {
    //		h->SetSubRange(i*n,n);
    //		b->SetSubRange(i*n,n);
    //		c->SetSubRange(i*n,n);

        Integrate(h(i*n,i*n+n-1),b(i*n,i*n+n-1), State2, 1/FS);

    //            Should be:  b = [0 , 1 , 3, 6, 10, 15, 21,... ]
    //            But becomes: b = [0, 0.1666, 1.1666, 3.166, 6.166, 10.166, 15.166, 21.166,... ]
    //            because of Simpson

        Differentiate(b(i*n,i*n+n-1),c(i*n,i*n+n-1),State,1/FS);

    //            Should be: c = [0,1,2,3,4,5,6....]
    //            But becomes: c = [0, 0.08333, 0.5833, 1.5, 2.5, 3.5, 4.5....]

      }
      c->SetFullRange(); //subranges were passed by reference!!
      DrawIt(c);
    }


    using Dew.Math;
    using Dew.Math.Units;
    using Dew.Signal;
    using Dew.Signal.Units;
    using Dew.Math.Tee;
    using Dew.Math.Editors;
    using Dew.Signal.Tee;

    private void button1_Click(object sender, EventArgs e)
    {
        int n,i;
        double FS = 1; //sampling frequency
        TIntegrateState State2 = new TIntegrateState();
        TDiffState State = new TDiffState();

        Vector h = MtxExpr.Ramp(30,0,1);
        Vector b = new Vector(h.Length);
        Vector c = new Vector(h.Length);

    //single block

        SignalUtils.Integrate(h, b, ref State2, 1.0/FS);
        SignalUtils.Differentiate(b, c, ref State, 1.0 / FS);
        TeeChart.DrawIt(c,"Processed in one block",false);

    //streaming

        //reset initial conditions
        State2 = new TIntegrateState();
        State = new TDiffState();

        n = h.Length / 10; //integer division (!)
        for (i = 0; i < 10; i++)
        {
            h.SetSubRange(i*n,n);
            b.SetSubRange(i*n,n);
            c.SetSubRange(i*n,n);
            SignalUtils.Integrate(h,b, ref State2, 1.0/FS);

    //            Should be:  b = [0 , 1 , 3, 6, 10, 15, 21,... ]
    //            But becomes: b = [0, 0.1666, 1.1666, 3.166, 6.166, 10.166, 15.166, 21.166,... ]
    //            because of Simpson

            SignalUtils.Differentiate(b,c, ref State,1.0/FS);

    //            Should be: c = [0,1,2,3,4,5,6....]
    //            But becomes: c = [0, 0.08333, 0.5833, 1.5, 2.5, 3.5, 4.5....]

        }
        c.SetFullRange();
        TeeChart.DrawIt(c,"Processed per partes, but same result",false);
    }