Note:  
The video below contains the playlist for all the videos in this series, which can be accessed via the playlist icon in the bottom-right corner of the embedded video frame once the video is playing. The first video in the series is loaded automatically

Note:  
This is an ongoing series of videos that will be updated every week. When a new video is posted we will update the publishing date of this article and the new video will be found at the end of the playlist

Intro to Matrices

Note:  
The video below contains the playlist for all the videos in this series, which can be accessed via the playlist icon in the bottom-right corner of the embedded video frame once the video is playing. The first video in the series is loaded automatically

Advanced Vectors

Note:  
The video below contains the playlist for all the videos in this series, which can be accessed via the playlist icon in the bottom-right corner of the embedded video frame once the video is playing. The first video in the series is loaded automatically

Intro to Vectors

Note:  
This article was originally published to GameDev.net back in 2002. It was revised by the original author in 2008 and published in the book Beginning Game Programming: A GameDev.net Collection, which is one of 4 books collecting both popular GameDev.net articls and new original content in print format.

Vectors

Vector Basics – What is a vector?

Vector Mathematics

Matrices

What is a Matrix anyway?

Matrix Arithmetic

Determinants

Matrix Transformations

Summary

References

Primitive description

Ray

• x represents any point on a line
• o represents some exact specified point on line (e.g. origin)
• d represents direction vector of line
• t is parameter within some specified range

class Ray
{
public:
    float4 origin;
    float4 direction;
	
    <<constructor>>
	
    <<methods>>
};

Plane

• n is unit vector in the plane normal direction
• x is every point on the plane
• p is our specified point on the plane

class Plane
{
public:
    float4 point;
    float4 normal;

    <<constructor>>
	
    <<methods>>
};

Sphere

• x, y and z are coordinates of point in sphere
• p represents sphere's center
• r represents sphere's radius

class Sphere
{
public:
    float4 center;
    float radius;
    
    <<constructor>>
    
    <<methods>>
};

Axis-aligned Box

• min_p stands for AABB minimum point
• max_p stands for AABB maximum point
• x represents points inside AABB

class AABB
{
public:
    float4 min_p;
    float4 max_p;
    
    <<constructor>>
    
    <<methods>>
};

Triangle

• x is every point in triangle
• A, B, C are points defining triangle
• a, ß, ? are so called barycentric coordinates (e.g. parameters in the equation above)

class Triangle
{
public:
    float4 p[3];
    
    <<constructor>>
        
    <<methods>>
};

Intersections

Ray-Plane intersection

bool IntrRayPlane(const Ray* r, const Plane* p, float* t)
{
    float denom = dot(p->normal, r->direction);
    
    // Test whether ray and plane aren't parallel
    if(fabsf(dotND) < EPSILON)
        return false;
    
    *t = -(dot(p->normal, r->origin - p->point) / denom);
    return ((*t) > 0.0f);
}

Ray-Sphere intersection

bool IntrRaySphere(const Ray* r, const Sphere* s, float* t)
{
    float rad2 = s->radius * s->radius;
    float4 L = s->center - r->origin;
    
    float tPX = dot(L, r->direction);
    if(tPX < 0.0)
        return false;
    
    float dsq = dot(L, L) - tPX * tPX;
    if(dsq > rad2)
        return false;
    
    float thit = sqrt(rad2 - dsq);
    *t = tPX - thit;
    if(*t < 0.0f)
        *t = tPX + thit;
    
    return ((*t) < 0.0f);
}    

Ray-AABB intersection

bool IntrRayAABB(const Ray* r, const AABB* b, float* enter, float* exit)
{
    float4 tmin = (b->minimum - r->origin) / r->direction;
    float4 tmax = (b->maximum - r->origin) / r->direction;
    
    float4 tnear = f4min(tmin, tmax);
    float4 tfar = f4min(tmin, tmax);
    
    *enter = max(max(tnear.x, 0.0f), max(tnear.y, tnear.z));
    *exit = min(tfar.x, min(tfar.y, tfar.z));
    
    return (*exit > 0.0f && *enter < *exit);
}

Ray-Triangle intersection

bool IntrRayTriangle( const Ray* r, const Triangle* t, float* thit, float* u, float* v)
{
    float4 e1 = t->v[1] - t->v[0];
    float4 e2 = t->v[2] - t->v[0];
    
    float4 pvec = cross(r->direction, e2);
    
    float det = dot(e1, pvec);
    if(det > -EPSILON && det < EPSILON)
        return false;
    
    float inv_det = 1.0f / det;
    
    float4 tvec = r->origin - t->v[0];
    *u = dot(tvec, pvec) * inv_det;
    if(*u < 0.0f || *u > 1.0f)
        return false;
    
    float4 qvec = cross(tvec, e1);
    *v = dot(r->direction, qvec) * inv_det;
    if(*v < 0.0f || *u + *v > 1.0f)
        return false;
    
    *thit = dot(e2, qvec) * inv_det;
    
    return (*thit > 0.0f);
}

Conclusion

Vector Math in Games

Concepts

Note:  In the world of programming programmers have utilized the vector type to represent both vectors in the traditional mathematical/physics sense as well as points or arbitrary n-tuplet units at their own discretion.  Be advised.

• Position - a vector represents an offset from your world origin point (0, 0, 0).
• Direction - vector looks very much like an arrow pointing at some direction. That is indeed what it can be used for. For example, if you have a vector that points south, you can make all your units go south.
• A special case of direction vector is a vector of length 1. It is called a normalized vector, or normal for short.
• A velocity vector can describe a movement. In that case, it is a difference in position over specific amount of time.

Remember the Basics - Vector Addition and Subtraction

Example: a Distance Between Objects

tree_position = (10, 10, 0)
my_position = (3, 3, 0)
# distance and direction you would need to move 
# to get exactly where the tree is
vector_to_tree = tree_position - my_position

Example: Velocity

position = (0, 10, 10) # position: 10 units in Y and Z direction
velocity = (500, 0, 0) # initial movement is 500 units in X direction over the next second

position += velocity # add velocity to position and update position

Concept: Simulation

position += velocity * delta

Example: Gravity, Air Resistance and Wind

# increase velocity of every object -2 down per second
gravity_vector = (0, -2, 0)

velocity += gravity_vector * delta # apply gravity effect
position += velocity * delta # update position

velocity += gravity_vector * delta # apply gravity effect
velocity *= 0.5 * delta # apply 0.5 slowdown per second
    
position += velocity * delta # update position

# modify kinetic energy / velocity
velocity += gravity_vector * delta # apply gravity effect
velocity *= 0.5 * delta # apply 0.5 slowdown per second
    
# add all forces
final_change_per_second = velocity + wind_force_per_second
    
# update position
position += final_change_per_second * delta

Concept: Direction

final_change = (B - A).normalize() * 3

Concept: Surface Plane

Concept: Dot Product

result = A dot B

• If dot product is 1, normals face the same direction.
• If dot product is 0, normals are perpendicular.
• If dot product is -1, normals face opposite directions.

angle = acos(A dot B)

Example: Light

• A normal vector which represent the direction of surface at this point.
• Position of the light.
• Position of this surface point.

distance_vec = light_pos - point_pos

light_direction = distance_vec.normalize()

brightness = surface_normal dot light_direction

varying vec3 surface_normal;
varying vec3 vertex_to_light_vector;

void main(void)
{
vec4 diffuse_color = vec3(1.0, 1.0, 1.0); // the color of surface - white
float diffuse_term = dot(surface_normal, normalize(vertex_to_light_vector));
gl_FragColor = diffuse_color * diffuse_term;
}

Example: Distance from a point to a plane

distance_to_a_plane = (point - plane_point) dot plane_normal;

Example: Is a point on a plane?

Example: Is a vector parallel to a plane?

Example: Line intersection with a plane

dot1 = SN dot (SP - P1)
dot2 = SN dot (P2 - P1)

u = (SN dot (SP - P1)) / (SN dot (P2 - P1))

• if u == 0, line is parallel the plane.
• if u <= 1 and u > 0, line intersects the plane.
• if u > 1, no intersection.

intersection point = (P2 - P1) * u

Concept: Cross Product

Example: Collision

temp_vector = change cross plane_normal

new_direction = temp_vector cross plane_normal

new_direction = (change cross plane_normal) cross plane_normal

Now what?

More in-depth reading

• Vector Math - Addition, subtraction, negation, normalization.
• Better Explained Dot Product - Follows similar style of explaining dot product by examples
• Dot Product and Cross Product - More in-depth explanation of dot and cross products
• What is a game loop? - Short explanation of the game loop
• More in-depth look at the different update loop types

struct Point {
  Vec2 Position;
  Vec2 OldPosition;
  Vec2 Acceleration;
};

class Physics {
  int PointCount;
  Point* Points[ MAX_VERTICES ];

  float Timestep;

public:
  void  UpdateVerlet();

  //Constructors, getters/setters etc. omitted
};

void Physics::UpdateVerlet() {
  for( int I = 0; I < PointCount; I++ ) {
    Point& P = *Points[ I ];

    Vec2 Temp = P.Position;
    P.Position += P.Position - P.OldPosition + P.Acceleration*Timestep*Timestep;
    P.OldPosition = Temp;
  }
}

struct Edge {
  Vertex* V1;
  Vertex* V2;

  float OriginalLength; //The length of the edge when it was created

  //Constructors etc. omitted
};

void Physics::UpdateEdges() {
  for( int I = 0; I < EdgeCount; I++ ) {
    Edge& E = *Edges[ I ];

    //Calculate the vector mentioned above
    Vec2 V1V2 = E.V2->Position - E.V1->Position; 

    //Calculate the current distance
    float V1V2Length = V1V2.Length(); 
    
    //Calculate the difference from the original length
    float Diff       = V1V2Length - E.OriginalLength; 
    
    V1V2.Normalize();

    //Push both vertices apart by half of the difference respectively 
    //so the distance between them equals the original length
    E.V1->Position += V1V2*Diff*0.5f; 
    E.V2->Position -= V1V2*Diff*0.5f;
  }
}

struct PhysicsBody {
  int VertexCount;
  int EdgeCount;

  Vertex* Vertices[ MAX_BODY_VERTICES ];
  Edge*   Edges   [ MAX_BODY_EDGES    ];

  void ProjectToAxis( Vec2& Axis, float& Min, float& Max );

  //Again, constructors etc. omitted
};

void PhysicsBody::ProjectToAxis( Vec2& Axis, float& Min, float& Max ) {
  float DotP = Axis*Vertices[ 0 ]->Position;

  //Set the minimum and maximum values to the projection of the first vertex
  Min = Max = DotP; 

  for( int I = 1; I < VertexCount; I++ ) {
    //Project the rest of the vertices onto the axis and extend 
    //the interval to the left/right if necessary
    DotP = Axis*Vertices[ I ]->Position; 

    Min = MIN( DotP, Min );
    Max = MAX( DotP, Max );
  }
}

bool Physics::DetectCollision( PhysicsBody* B1, PhysicsBody* B2 ) {
  //Just a fancy way of iterating through all of the edges of both bodies at once
  for( int I = 0; I < B1->EdgeCount + B2->EdgeCount; I++ ) { 
    Edge* E;

    if( I < B1->EdgeCount )
      E = B1->Edges[ I ];
    else
      E = B2->Edges[ I - B1->EdgeCount ];
    
    //Calculate the axis perpendicular to this edge and normalize it
    Vec2 Axis( E->V1->Position.Y - E->V2->Position.Y, E->V2->Position.X - E->V1->Position.X ); 
    Axis.Normalize();

    float MinA, MinB, MaxA, MaxB; //Project both bodies onto the perpendicular axis
    B1->ProjectToAxis( Axis, MinA, MaxA );
    B2->ProjectToAxis( Axis, MinB, MaxB );

    //Calculate the distance between the two intervals - see below
    float Distance = IntervalDistance( MinA, MaxA, MinB, MaxB ); 

    if( Distance > 0.0f ) //If the intervals don't overlap, return, since there is no collision
      return false;
  }

  return true; //There is no separating axis. Report a collision!
}

float Physics::IntervalDistance( float MinA, float MaxA, float MinB, float MaxB ) {
  if( MinA < MinB )
    return MinB - MaxA;
  else
    return MinA - MaxB;
}

class Physics {
  struct {
    float Depth;
    Vec2  Normal;
  } CollisionInfo;

  //Everything else omitted
}

bool Physics::DetectCollision( PhysicsBody* B1, PhysicsBody* B2 ) {
  float MinLength = 10000.0f; //Initialize the length of the collision vector to a relatively large value
  for( int I = 0; I < B1->EdgeCount + B2->EdgeCount; I++ ) {
    Edge* E;

    if( I < B1->EdgeCount )
      E = B1->Edges[ I ];
    else
      E = B2->Edges[ I - B1->EdgeCount ];

    Vec2 Axis( E->V1->Position.Y - E->V2->Position.Y, E->V2->Position.X - E->V1->Position.X );
    Axis.Normalize();

    float MinA, MinB, MaxA, MaxB;
    B1->ProjectToAxis( Axis, MinA, MaxA );
    B2->ProjectToAxis( Axis, MinB, MaxB );

    float Distance = IntervalDistance( MinA, MaxA, MinB, MaxB );

    if( Distance > 0.0f )
      return false;
      
    //If the intervals overlap, check, whether the vector length on this 
    //edge is smaller than the smallest length that has been reported so far
    else if( abs( Distance ) < MinDistance ) { 
      MinDistance = abs( Distance );

      CollisionInfo.Normal = Axis; //Save collision information for later
    }
  }

  CollisionInfo.Depth = MinDistance;

  return true; //There is no separating axis. Report a collision!
}

struct {
  float Depth;
  Vec2  Normal;

  Edge*   E;
  Vertex* V;
} CollisionInfo;

bool Physics::DetectCollision( PhysicsBody* B1, PhysicsBody* B2 ) {
  float MinDistance = 10000.0f;
  for( int I = 0; I < B1->EdgeCount + B2->EdgeCount; I++ ) { //Same old
    Edge* E;

    if( I < B1->EdgeCount )
      E = B1->Edges[ I ];
    else
      E = B2->Edges[ I - B1->EdgeCount ];

    Vec2 Axis( E->V1->Position.Y - E->V2->Position.Y, E->V2->Position.X - E->V1->Position.X );
    Axis.Normalize();

    float MinA, MinB, MaxA, MaxB;
    B1->ProjectToAxis( Axis, MinA, MaxA );
    B2->ProjectToAxis( Axis, MinB, MaxB );

    float Distance = IntervalDistance( MinA, MaxA, MinB, MaxB );

    if( Distance > 0.0f )
      return false;
    else if( abs( Distance ) < MinDistance ) {
      MinDistance = abs( Distance );

      CollisionInfo.Normal = Axis;
      CollisionInfo.E      = E; //Store the edge, as it is the collision edge
    }
  }

  CollisionInfo.Depth = MinDistance;

  //Ensure that the body containing the collision edge lies in 
  //B2 and the one containing the collision vertex in B1
  if( CollisionInfo.E->Parent != B2 ) { 
    PhysicsBody* Temp = B2;
    B2 = B1;
    B1 = Temp;
  }

  //This is needed to make sure that the collision normal is pointing at B1
  int Sign = SGN( CollisionInfo.Normal*( B1->Center - B2->Center ) ); 

  //Remember that the line equation is N*( R - R0 ). We choose B2->Center 
  //as R0; the normal N is given by the collision normal

  if( Sign != 1 )
    CollisionInfo.Normal = -CollisionInfo.Normal; //Revert the collision normal if it points away from B1


  float SmallestD = 10000.0f; //Initialize the smallest distance to a high value
  for( int I = 0; I < B1->VertexCount; I++ ) {
    //Measure the distance of the vertex from the line using the line equation
    float Distance = CollisionInfo.Normal*( B1->Vertices[ I ]->Position - B2->Center ); 

    //If the measured distance is smaller than the smallest distance reported 
    //so far, set the smallest distance and the collision vertex
    if( Distance < SmallestD ) { 
      SmallestD = Distance;
      CollisionInfo.V = B1->Vertices[ I ];
    }
  }

  return true;
}

void Physics::ProcessCollision() {
  Vec2 CollisionVector = CollisionInfo.Normal*CollisionInfo.Depth;

  CollisionInfo.V->Position += CollisionVector*0.5f;
}

Vertex* E1 = CollisionInfo.E->V1;
Vertex* E2 = CollisionInfo.E->V2;

float T = ( CollisionInfo.V->Position.X - CollisionVector.X - E1->Position.X )/(  E2->Position.X - E1->Position.X );

float T;
if( abs( E1->Position.X - E2->Position.X ) > abs( E1->Position.Y - E2->Position.Y ) )
  T = ( CollisionInfo.V->Position.X - CollisionVector.X - E1->Position.X )/(  E2->Position.X - E1->Position.X );
else
  T = ( CollisionInfo.V->Position.Y - CollisionVector.Y - E1->Position.Y )/(  E2->Position.Y - E1->Position.Y );

void Physics::ProcessCollision() {
  Vec2 CollisionVector = CollisionInfo.Normal*CollisionInfo.Depth;

  Vertex* E1 = CollisionInfo.E->V1;
  Vertex* E2 = CollisionInfo.E->V2;

  float T;
  if( abs( E1->Position.X - E2->Position.X ) > abs( E1->Position.Y - E2->Position.Y ) )
    T = ( CollisionInfo.V->Position.X - CollisionVector.X - E1->Position.X )/(  E2->Position.X - E1->Position.X);
  else
    T = ( CollisionInfo.V->Position.Y - CollisionVector.Y - E1->Position.Y )/(  E2->Position.Y - E1->Position.Y);

  float Lambda = 1.0f/( T*T + ( 1 - T )*( 1 - T ) );

  E1->Position -= CollisionVector*( 1 - T )*0.5f*Lambda;
  E2->Position -= CollisionVector*      T  *0.5f*Lambda;

  CollisionInfo.V->Position += CollisionVector*0.5f;
}

void Physics::Update() {
  UpdateForces();
  UpdateVerlet();
  IterateCollisions();
}

Introduction

Someone asked me in the #gamedev IRC channel about how to make a 2d vehicle simulator. Instead of spending all day trying to explain the concepts to him, I decided just to write this tutorial.
Please bear with me, this is my first tutorial.

So as I mentioned we’re going to be learning how to make a basic 2d vehicle simulator. We’re going to do it in C# and try do use as few hacks as possible. I’ve broken the process
down into three steps. First, we will learn how to setup a basic game application in C#.NET and how to draw some basic graphics (emphasis on basic.) Next, we will learn how to create a rigid body
simulator using a simple Euler integrator with a variable time step. And last but not least, we will calculate vehicle forces simulating the tire patch contacting the road. And that’s all there
is to it! Let’s get started.

Math Requirements

Phase One

Phase one, as I mentioned, is to create the renderer; something graphical so we can actually see what our simulation is doing. This will make it a lot easier to debug. Create a windows form
project in C# and place a picturebox control on it (name it "screen"). This control is where we will display our simulation. We could just start drawing to this screen but we’re going to be
using double buffering as well to avoid flicker, so we need to create the back buffer now. That bit of code looks like this.

Graphics graphics; //gdi+

Bitmap backbuffer;

Size buffersize;



//intialize rendering

private void Init(Size size)

{

  //setup rendering device

  buffersize = size;

  backbuffer = new Bitmap(buffersize.Width, buffersize.Height);

  graphics = Graphics.FromImage(backbuffer);

}

//main rendering function

private void Render(Graphics g)

{

  //clear back buffer

  graphics.Clear(Color.Black);



  //draw to back buffer

  graphics.DrawLine(new Pen(Color.Yellow), 1, 0, 1, 5);



  //present back buffer

  g.DrawImage(backbuffer, new Rectangle(0, 0, buffersize.Width, buffersize.Height), 0, 0, buffersize.Width, buffersize.Height, GraphicsUnit.Pixel);

}

Now by default, the graphics of a picturebox control has the origin in the top-left corner, and extends downward for +y and to the right for +x. This is highly unnatural for most cases. In
addition to that, it has extremely large units. Since we will be simulating in the metric system, I recommend introducing a scale factor to scale up the simulation and make it much more visible. The
transformation looks like this and takes place after Graphics.Clear() is called.

graphics.ResetTransform();

graphics.ScaleTransform(screenScale, -screenScale);

graphics.TranslateTransform(buffersize.Width / 2.0f / screenScale, -buffersize.Height / 2.0f / screenScale);

Now the line should draw starting right at the center of the screen.

Forms Wiring

Up until now, I havn’t explained how to connect all the functions. The first thing you’ll need to do is call the Render function from your on_paint event.
Next, you’ll need to create a function that gets called continously to update the simulation. It is preferred to call this function on the Application_Idle event. So create an
event handler for Application_Idle and have it call your DoFrame function. Inside this function you’ll need to

1. Process input
2. Update the simulation
3. Invalidate the screen control

The Timer

Since we don’t know how often our DoFrame function will be getting called, we need to code everything to handle a variable time step. To utilize this we must measure the time
between DoFrame calls. So I’ll introduce the timer which, very simply, queries the number of milliseconds that have passed since the computer was turned on. So we store this number
every frame and on a subsequent frame we compute the difference, which gives us the amount of time that has passed since the last frame. Here is my very simple timer object. Note: you will
need to call GetETime in your intialize function in order to clear the timer, otherwise the first call to it will return the amount of time that has passed since the computer was turned
on.

class Timer

{

  //store last time sample

  private int lastTime = Environment.TickCount;

  private float etime;



  //calculate and return elapsed time since last call

  public float GetETime()

  {

    etime = (Environment.TickCount - lastTime) / 1000.0f;

    lastTime = Environment.TickCount;



    return etime;

  }

}

Conclusion of Phase One

So up until now we’ve covered: setting up a rendering surface using GDI, wiring a form to process a game loop and draw it to the screen, and computing the time that has passed since the last
frame. Our application looks like this:

using System;

using System.Collections.Generic;

using System.ComponentModel;

using System.Data;

using System.Drawing;

using System.Drawing.Drawing2D;

using System.Text;

using System.Windows.Forms;



namespace racing_simulation_2d

{

  //our main application form

  public partial class frmMain : Form

  {

    //graphics

    Graphics graphics; //gdi+

    Bitmap backbuffer;

    Size buffersize;

    const float screenScale = 3.0f;

    Timer timer = new Timer();



    //keyboard controls

    bool leftHeld = false, rightHeld = false;

    bool upHeld = false, downHeld = false;



    //vehicle controls

    float steering = 0; //-1 is left, 0 is center, 1 is right

    float throttle = 0; //0 is coasting, 1 is full throttle

    float brakes = 0; //0 is no brakes, 1 is full brakes



    public frmMain()

    {

      InitializeComponent();

      Application.Idle += new EventHandler(ApplicationIdle);



      screen.Paint += new PaintEventHandler(screen_Paint);

      this.KeyDown += new KeyEventHandler(onKeyDown);

      this.KeyUp += new KeyEventHandler(onKeyUp);



      Init(screen.Size);

    }



    //intialize rendering

    private void Init(Size size)

    {

      //setup rendering device

      buffersize = size;

      backbuffer = new Bitmap(buffersize.Width, buffersize.Height);

      graphics = Graphics.FromImage(backbuffer);



      timer.GetETime(); //reset timer

    }



    //main rendering function

    private void Render(Graphics g)

    {

      //clear back buffer

      graphics.Clear(Color.Black);

      graphics.ResetTransform();

      graphics.ScaleTransform(screenScale, -screenScale);

      graphics.TranslateTransform(buffersize.Width / 2.0f /

          screenScale, -buffersize.Height / 2.0f / screenScale);



      //draw to back buffer

      DrawScreen();



      //present back buffer

      g.DrawImage(backbuffer, new Rectangle(0, 0, buffersize.Width,

       buffersize.Height), 0, 0, buffersize.Width, 

       buffersize.Height, GraphicsUnit.Pixel);

    }



    //draw the screen

    private void DrawScreen()

    {

      //draw our simulation here

    }



    //process game logic

    private void DoFrame()

    {

      //get elapsed time since last frame

      float etime = timer.GetETime();



      //process input

      ProcessInput();



      ////////////////////////////////

      //integrate our simulation here

      ////////////////////////////////



      //redraw our screen

      screen.Invalidate();

    }





    //process keyboard input

    private void ProcessInput()

    {

      if (leftHeld)

        steering = -1;

      else if (rightHeld)

        steering = 1;

      else

        steering = 0;



      if (upHeld)

        throttle = 1;

      else

        throttle = 0;



      if (downHeld)

        brakes = 1;

      else

        brakes = 0;

    }



    private void onKeyDown(object sender, KeyEventArgs e)

    {

      switch (e.KeyCode)

      {

        case Keys.Left:

          leftHeld = true;

          break;

        case Keys.Right:

          rightHeld = true;

          break;

        case Keys.Up:

          upHeld = true;

          break;

        case Keys.Down:

          downHeld = true;

          break;

        default: //no match found

          return; //return so handled dosnt get set

      }



      //match found

      e.Handled = true;

    }



    private void onKeyUp(object sender, KeyEventArgs e)

    {

      switch (e.KeyCode)

      {

        case Keys.Left:

          leftHeld = false;

          break;

        case Keys.Right:

          rightHeld = false;

          break;

        case Keys.Up:

          upHeld = false;

          break;

        case Keys.Down:

          downHeld = false;

          break;

        default: //no match found

          return; //return so handled dosnt get set

      }



      //match found

      e.Handled = true;

    }



    //rendering - only when screen is invalidated

    private void screen_Paint(object sender, PaintEventArgs e)

    {

      Render(e.Graphics);

    }



    //when the os gives us time, run the game

    private void ApplicationIdle(object sender, EventArgs e)

    {

      // While the application is still idle, run frame routine.

      DoFrame();

    }



    private void MenuExit_Click(object sender, EventArgs e)

    {

      this.Close();

    }

  }



  //keep track of time between frames

  class Timer

  {

    //store last time sample

    private int lastTime = Environment.TickCount;

    private float etime;



    //calculate and return elapsed time since last call

    public float GetETime()

    {

      etime = (Environment.TickCount - lastTime) / 1000.0f;

      lastTime = Environment.TickCount;



      return etime;

    }

  }

}

Phase Two - Rigid Body Simulation

Ok, now we’re getting into some good stuff here. Let's put everything we just covered on the back burner now and talk about some physics. We’re going to be using a very simple Euler
integration method. Basically, each frame we accumulate a bunch of forces (in our case from each wheel of the vehicle) and calculate the resultant acceleration, which is in the form of A=F/M (the
same as F=MA, Newton’s second law of motion). We use this to modify Newton’s first law of motion, “an object in motion stays in motion…” So we calculate our A, and we
integrate it into our V. Without an A, V would be constant, hence staying in motion, if no forces should act on it. Newton's third law gets applied in the form that any potential force the vehicle is
applying to the ground, gets applied in the opposite direction to the vehicle (I'll explain this in the vehicle section). This topic is much easier to explain with symbols. So, P is our vehicle
position, V is its linear velocity, F is the net force acting on it, M is its mass, A is the resultant acceleration, and T is the time step (the value our timer gave us from the last frame).

  A = F / M

  V = V + A * T

  P = P + V * T

This is a basic linear rigid body simulator. Each frame, we total up some F, integrate it, and then zero out F to restart the accumulation the next frame. Now let’s talk about rotation. The
angular case is nearly identical to the linear case (especially in 2D). Instead of P we have an Angle, instead of V we have an Angular Velocity, instead of F we have a torque, and instead of M we
have inertia. So the angular model looks like this

  AngA = Torque / Inertia

  AngV = AngV + AngA * T

  Angle = Angle + AngV * T

This is what that bit of code looks like. % is the cross product operator for my vector class.

public void AddForce(Vector worldForce, Vector worldOffset)

{

  //add linar force

  m_forces += worldForce;

  //and it's associated torque

  m_torque += worldOffset % worldForce;

}

Code Dump

Here is my rigid body object. You’ll notice all the properties I mentioned above. It has a Draw function which will draw its rectangle to the provided graphics object. It has an
AddForce function, a space conversion method, to and from world space (very handy), and a function that returns the velocity of a point on the body (in world space). This point velocity
is a combination of the linear velocity and the angular velocity. But the angular velocity is multiplied by the distance the point is from the center of rotation and perpendicular to its offset
direction. So to kill two birds with one stone, I simply find the orthogonal vector to the point offset and multiply it by the angular velocity (then add the linear velocity.)

One thing you may be curious about is how I calculate the inertia value. That is a generalized formula I found at this
link.

//our simulation object

class RigidBody

{

  //linear properties

  private Vector m_position = new Vector();

  private Vector m_velocity = new Vector();

  private Vector m_forces = new Vector();

  private float m_mass;



  //angular properties

  private float m_angle;

  private float m_angularVelocity;

  private float m_torque;

  private float m_inertia;



  //graphical properties

  private Vector m_halfSize = new Vector();

  Rectangle rect = new Rectangle();

  private Color m_color;



  public RigidBody()

  { 

    //set these defaults so we dont get divide by zeros

    m_mass = 1.0f; 

    m_inertia = 1.0f; 

  }



  //intialize out parameters

  public void Setup(Vector halfSize, float mass, Color color)

  {

    //store physical parameters

    m_halfSize = halfSize;

    m_mass = mass;

    m_color = color;

    m_inertia = (1.0f / 12.0f) * (halfSize.X * halfSize.X)

       * (halfSize.Y * halfSize.Y) * mass;



    //generate our viewable rectangle

    rect.X = (int)-m_halfSize.X;

    rect.Y = (int)-m_halfSize.Y;

    rect.Width = (int)(m_halfSize.X * 2.0f);

    rect.Height = (int)(m_halfSize.Y * 2.0f);

  }



  public void SetLocation(Vector position, float angle)

  {

    m_position = position;

    m_angle = angle;

  }



  public Vector GetPosition()

  {

    return m_position;

  }



  public void Update(float timeStep)

  {

    //integrate physics

    //linear

    Vector acceleration = m_forces / m_mass;

    m_velocity += acceleration * timeStep;

    m_position += m_velocity * timeStep;

    m_forces = new Vector(0,0); //clear forces



    //angular

    float angAcc = m_torque / m_inertia;

    m_angularVelocity += angAcc * timeStep;

    m_angle += m_angularVelocity * timeStep;

    m_torque = 0; //clear torque

  }



  public void Draw(Graphics graphics, Size buffersize)

  {

    //store transform, (like opengl's glPushMatrix())

    Matrix mat1 = graphics.Transform;



    //transform into position

    graphics.TranslateTransform(m_position.X, m_position.Y);

    graphics.RotateTransform(m_angle/(float)Math.PI * 180.0f);



    try

    {

      //draw body

      graphics.DrawRectangle(new Pen(m_color), rect);



      //draw line in the "forward direction"

      graphics.DrawLine(new Pen(Color.Yellow), 1, 0, 1, 5);

    }

    catch(OverflowException exc)

    {

      //physics overflow :(

    }  



    //restore transform

    graphics.Transform = mat1;

  }



  //take a relative vector and make it a world vector

  public Vector RelativeToWorld(Vector relative)

  {

    Matrix mat = new Matrix();

    PointF[] vectors = new PointF[1];



    vectors[0].X = relative.X;

    vectors[0].Y = relative.Y;



    mat.Rotate(m_angle / (float)Math.PI * 180.0f);

    mat.TransformVectors(vectors);



    return new Vector(vectors[0].X, vectors[0].Y);

  }



  //take a world vector and make it a relative vector

  public Vector WorldToRelative(Vector world)

  {

    Matrix mat = new Matrix();

    PointF[] vectors = new PointF[1];



    vectors[0].X = world.X;

    vectors[0].Y = world.Y;



    mat.Rotate(-m_angle / (float)Math.PI * 180.0f);

    mat.TransformVectors(vectors);



    return new Vector(vectors[0].X, vectors[0].Y);

  }



  //velocity of a point on body

  public Vector PointVel(Vector worldOffset)

  {

    Vector tangent = new Vector(-worldOffset.Y, worldOffset.X);

    return tangent * m_angularVelocity + m_velocity;

  }



  public void AddForce(Vector worldForce, Vector worldOffset)

  {

    //add linar force

    m_forces += worldForce;

    //and it's associated torque

    m_torque += worldOffset % worldForce;

  }

}

Testing

To make sure your rigid body works, instantiate one in your Init() function and apply a force with some offset in the DoFrame function. If you apply a constant
worldOffset, the body will continue to accelerate its angular velocity. If you take your offset and run it through the RelativeToWorld function, the body will angularly accelerate in one
direction and then come back the other way, like a pendulum as the point the force is applied to changes. Play around with this for a while, this has to work and make sense in order for the next
section to work.

Phase Three - The Vehicle

Assuming everything has gone well above, you should have a rigid body actor in your scene that you can apply forces to and watch move around. Now all that’s left is to calculate these forces
in a way that will simulate a vehicle. For that we are going to need a vehicle object. I recommend deriving directly from you rigid body object since the chassis is essentially a rigid body. In
addition to that we will need to construct a “wheel" object. This wheel will handle the steering direction of each wheel, the velocity the wheel is spinning, and calculate the forces that that
particular wheel applies to the chassis (all in vehicle space). Since our wheel is known to be constrained to the vehicle, we don’t need to simulate it as another rigid body (though you could,
but not in the 2D case.) We will simply duplicate the angular properties of the rigid body in the wheel object.

So we’ll need: Wheel Velocity, Wheel Inertia, and and Wheel Torque. We’ll also need the relative offset of the wheel in the vehicle space, and the angle the wheel is facing (this is
constant for the back wheels, unless you want 4 wheel steering.) Just like the rigid body, the wheel's torque function acts as an accumulator, we add torques to it and after it gets integrated the
torque is zeroed out. The AddTorque function is where you will apply a wheel torque from either the transmission (to make you go) or from the brakes (to make you stop). Internally the
wheel will generate a torque caused by the friction on the road.

The wheel object also needs a SetSteering function. This function calculates two vectors: an effective Side Direction, and an effective Forward Direction (both in vehicle space) that
the tire patch will act on. The force applied on the tire by the ground acting in the side direction will directly translate into the chassis. Meanwhile the force acting in the forward direction will
not only act on the chassis, but it will induce a rotation of the tire. Here is the SetSteering function; you will see I used the Drawing2D.Matrix to transform the initial
forward and side vectors by the steering angle (I had to convert the vectors to “points" in order to transform them by the matrix.)

public void SetSteeringAngle(float newAngle)

{

  Matrix mat = new Matrix();

  PointF[] vectors = new PointF[2];



  //foward vector

  vectors[0].X = 0;

  vectors[0].Y = 1;

  //side vector

  vectors[1].X = -1;

  vectors[1].Y = 0;



  mat.Rotate(newAngle / (float)Math.PI * 180.0f);

  mat.TransformVectors(vectors);



  m_forwardAxis = new Vector(vectors[0].X, vectors[0].Y);

  m_sideAxis = new Vector(vectors[1].X, vectors[1].Y);

}

Force Calculation

So, if the vehicle is sitting there not moving with its front wheels turned, and you push it, a force will be generated in the opposite direction you push. This force gets projected onto these two
directions. If the wheels were straight there would be no side force. So the vehicle would simply roll forward. But since the wheels are turned, there is a bit of the force that acts in the
“effective side direction” so we apply an opposite force to the chassis. This is what causes you to turn when you steer the wheels. To get this force that gets projected onto the two
directions, we need to first determine the velocity difference between the tire patch and the road. If the wheel is spinning at the same speed the ground is wizzing by, then there is effectively no
force acting on the vehicle. But as soon as you slam on the brakes and stop the wheel, there is a huge velocity difference and this is what causes the force that stops your car.

So here is the process broken down into 6 steps, for each wheel.

Step 1, calculate the effective direction vectors (with steering function).

Step 2, calculate velocity difference. The ground speed is determined via the “PointVel” function on the rigidbody, given the current wheel’s world offset.

Step 3, project this velocity onto the two effective directions.

Step 4, generate an equal and opposite force for the two direction and call this the “response force”. This is what gets added to the chassis for each wheel.

Step 5, calculate the torque that the forward response force created on the wheel, and add this to the wheel torque.

Step 6, integrate the wheel torques into the wheel velocity.

That bit of code looks like this:

public Vector CalculateForce(Vector relativeGroundSpeed, float timeStep)

{

  //calculate speed of tire patch at ground

  Vector patchSpeed = -m_forwardAxis * m_wheelSpeed *

    m_wheelRadius;



  //get velocity difference between ground and patch

  Vector velDifference = relativeGroundSpeed + patchSpeed;



  //project ground speed onto side axis

  float forwardMag = 0;

  Vector sideVel = velDifference.Project(m_sideAxis);

  Vector forwardVel = velDifference.Project(m_forwardAxis, out forwardMag);



  //calculate super fake friction forces

  //calculate response force

  Vector responseForce = -sideVel * 2.0f;

  responseForce -= forwardVel;



  //calculate torque on wheel

  m_wheelTorque += forwardMag * m_wheelRadius;



  //integrate total torque into wheel

  m_wheelSpeed += m_wheelTorque / m_wheelInertia * timeStep;



  //clear our transmission torque accumulator

  m_wheelTorque = 0;



  //return force acting on body

  return responseForce;

}

Almost Done!

We’re in the home stretch here now. Now we have a way to calculate the force each wheel generates on the chassis. Every frame, all we have to do is set our transmission and brake torques,
our steering angle, calculate each wheel force, add these to the chassis, and integrate the rigid body. Badaboom badabing, vehicle done! :)

Conclusion

Here is the entire source code for the project. If you have any questions or comments please
feel free to post them here and either I can make things more clear or maybe someone else could offer some better expertise. If you’d like you can email me at Kincaid05 on google's fine
emailing service.

Thanks for reading and I hope this was informative.
-Matt Kincaid

• The objects that the algorithm is dealing with are convex³ 3D objects only. It is true that the algorithm will work for some concave⁴ objects but care must be taken⁵.
• Detection of any pair of faces colliding is a well-defined and solved problem.
• Both objects have a pre-calculated center point.
• The accuracy constant k is set to be (min/max estimates)
• Object A has exactly n faces while object B has m faces.

• Optimization: check collision between object A's bounding sphere and object B's bounding sphere.
• Find the closest k faces of object A relative to object B.
• Calculate the geometric center of the new selection and the bounding sphere radius.
• Find the closest k faces of object B relative to object A's new selection of k faces.
• Calculate the geometric center of object B's new selection of faces and their maximal bounding sphere radius.
• Optimization: check collision between spheres to determine if there is even a chance for face collisions.
• Sort the two sets of faces by increasing distance (optional, might be replaces by a good insertion algorithm).
• Test the two sets of faces against each other, starting with the closest pairs of faces.

• Object B can be considered as a point object located at its center (object centers have been shown in step 1).

• Using both object's positions, calculate the relative position vector. Assuming is the center of object A (x,y,z - a 3D vector) and the center of object B, the vector connecting the center points is then:
• Start looping through all faces of object A, find each face center relative to object A's center. Total time complexity is O(n). It takes less operations to calculate the faces centers than to do it in the beginning and apply transformations on them.
• For each face, find the vector connecting its center with object B's center, using the simple relation: , where is the vector connecting the center of face i with the center of object B, the relative positions of the two object's centers and is the position of the face relative to object A's center (the object it belongs to). Of course, i = 0,1,2,3...,n.
• For each , find the vector's size: and store into selection*. Of course, the size itself is irrelevant, as is the direction of the vector. The only important quantity is the vector's size squared. Note that there is no need to take the square root because if a certain face's squared vector size is larger than another face's squared vector size, the same holds for the actual sizes themselves.
• Depending on the insertion technique used to insert the distances (squared sizes of the relative position vectors), find the smallest k distances which should be straightforward if the insertion was done in an efficient way, and remember their corresponding faces.

• Calculate the relative position vector between object B's center and the center of the new selection of k faces found earlier. Assuming is the center of object B and the center of the selection in object A, the vector connecting the center points is then:
• Start looping through all faces of object B, find each face center relative to the center of the selection in object A. Total time complexity is O(m). It takes less operations to calculate the faces centers than to do it in the beginning and apply transformations on them.
• For each face, find the vector connecting its center with the center of the selection in object A, using the simple relation: , where is the vector connecting the center of face i with the center of the selection in object A, the relative positions of the center of object B and the center of the selection in object A and is the position of the face relative to object B's center (the object it belongs to). Of course, i = 0,1,2,3...,m.
• For each , find the vector's size: and store into selection*.
• Depending on the insertion technique used to insert the distances (squared sizes of the relative position vectors), find the smallest k distances which should be straightforward if the insertion was done in an efficient way, and remember their corresponding faces.

• The algorithm supports detection of mostly convex objects. Some concave objects will work as well but there are bound to be degenerate cases that cause the algorithm to fail. It hasn't been proven mathematically that there are any cases that will cause failure but assuming there aren't is not a good idea.
• Moving objects in computer simulations are by nature problematic because time in computer systems is discrete. In numerical approximations, time is discrete and therefore object positions are calculated in time steps and are really "teleporting" in small steps to create the illusion of motion. The problem arises when the velocity of objects is very large and/or objects are very small. It is quite possible that an object will be close to another object while not colliding, yet at the next time step will be half intersecting with the other object. The problem is then what faces to test. After all, the "closest" faces the object previously had are now inside the other close object. Even more so, it's possible that if the velocity is large enough, the object might pass right through the other object without any way of us detecting this.
• In a similar fashion, rotating objects pose another problem. If objects have an angular velocity (or momentum, whichever you prefer) it's possible that before the collision some faces are the closest selection while after a small time step, other faces should be the closest selection. For example, consider a normal cube floating above a table at a small height. The closest faces to the table are naturally the two faces (triangles) that comprise the base of the cube. If the cube is about to rotate 45 degrees in one time step (very fast rotation) it wouldn't be correct checking the base of the cube against the table.

• Find the closest object in collection A relative to the entire collection B. This implies each object in A has to have a center (this is mandatory for OFG anyway) and the entire collection has to have a center as well (averaging the centers of the objects, should be a part of initialization).
• Find the closest object in collection B relative to the object found earlier in collection A.
• Feed the two objects found to the OFG algorithm.

• Calculate the translation of object A after a time unit elapses.
• Feed the two objects into the OFG algorithm, forgoing step 6 completely and without the bounding sphere calculations in steps 3 and 5. Also, there is no need to perform step 8 (the final step) of the OFG algorithm just yet.
• Connect each vertex in selection A with the same vertex on the virtual selection of A that is translated by the amount calculated in step 1. This gives a vector that represents the edge to test against object B.
• Test the newly created edges against the selection in object B found by the OFG algorithm. They are the most likely to be intersecting object B.
• Test the original faces in selection A against the selection in B, as was intended in step 8 of the OFG algorithm.
• Test the translated virtual selection of A against the selection of B using step 8 of the OFG algorithm.

• Using the translation found in step 1 earlier and given that it's easily possible to calculate the relative velocity of the two objects, it's possible to estimate the time of collision (therefore obtaining the delta time needed for the objects to collide).
• Moving half of the distance (half a time unit) and checking for collisions. If collision occurred, great. If not, move another half, so on and so forth. This has a time complexity of log(n) times of performing all the collision tests.

• Solve the plane equation of the second face for each vertex in the first face. If all vertices produce the same sign, the first face definitely doesn't pierce the plane in question. Creating the plane equation is easy given two vectors on that plane (three vertices of one face can give two vectors). If any vertex gives the other sign, there is a possibility of collision.

• Calculate the intersection points of the edges found piercing the plane in stage one using the second face's plane equation.
• For each edge in the second face, take the cross product of the edge vector with the normal to the face, giving a perpendicular plane.
• The points must lie either on the plane or on either side. Using the same logic as in stage one can determine this for each edge plane.
• If the point lies in the same part of space ("side") for all planes corresponding to the edges of the face, the point is within the face and therefore a collision has occurred
• If both points lie outside the planes of the edges (not all signs are the same for each point checked against the edge planes), there is definitely no collision.

So, what is a Quaternion?

Q = xi + yj + zk + w  

Q = [ w, [b]v[/b] ]

struct quaternion
{
  double x, y, z, w;
};

Normalizing a quaternion

|Q| = sqrt( w^2 + x^2 + y^2 + z^2)

double length(quaternion quat)
{
  return sqrt(quat.x * quat.x + quat.y * quat.y +
              quat.z * quat.z + quat.w * quat.w);
}

Q* = Q/|Q| = [w/|Q|, [b]v[/b]/|Q|]

quaternion normalize(quaternion quat)
{
  double L = length(quat);

  quat.x /= L;
  quat.y /= L;
  quat.z /= L;
  quat.w /= L;

  return quat;
}

The conjugate of a quaternion

Q' = [ w, -[b]v[/b] ]

quaternion conjugate(quaternion quat)
{
  quat.x = -quat.x;
  quat.y = -quat.y;
  quat.z = -quat.z;
  return quat;
}

Multiplying quaternions

C = A * B

C.x = | A.w*B.x + A.x*B.w + A.y*B.z - A.z*B.y |
C.y = | A.w*B.y - A.x*B.z + A.y*B.w + A.z*B.x |
C.z = | A.w*B.z + A.x*B.y - A.y*B.x + A.z*B.w |
C.w = | A.w*B.w - A.x*B.x - A.y*B.y - A.z*B.z |

quaternion mult(quaternion A, quaternion B)
{
  quaternion C;

  C.x = A.w*B.x + A.x*B.w + A.y*B.z - A.z*B.y;
  C.y = A.w*B.y - A.x*B.z + A.y*B.w + A.z*B.x;
  C.z = A.w*B.z + A.x*B.y - A.y*B.x + A.z*B.w;
  C.w = A.w*B.w - A.x*B.x - A.y*B.y - A.z*B.z;
  return C;
}

The quaternion camera

V = [0, View]

vector A = [x, y, z] 
R.x = A.x * sin(theta/2)
R.y = A.y * sin(theta/2)
R.z = A.z * sin(theta/2)
R.w = cos(theta/2)

W = R * V * R'

NewView = [W.x W.y W.z]

void Camera::SetViewByMouse(void)
{
  // the coordinates of our mouse coordinates
  int MouseX, MouseY;

  // the middle of the screen in the x direction
  int MiddleX = SCREENWIDTH/2;

  // the middle of the screen in the y direction
  int MiddleY = SCREENHEIGHT/2;

  // vector that describes mouseposition - center
  Vector MouseDirection(0, 0, 0);

  // static variable to store the rotation about the x-axis, since
  // we want to limit how far up or down we can look.
  // We don't need to cap the rotation about the y-axis as we
  // want to be able to turn around 360 degrees
  static double CurrentRotationAboutX = 0.0;

  // The maximum angle we can look up or down, in radians
  double maxAngle = 1;

  // This function gets the position of the mouse
  SDL_GetMouseState(&MouseX, &MouseY);

  // if the mouse hasn't moved, return without doing
  // anything to our view
  if((MouseX == MiddleX) && (MouseY == MiddleY))
    return;

  // otherwise move the mouse back to the middle of the screen
  SDL_WarpMouse(MiddleX, MiddleY);

  // get the distance and direction the mouse moved in x (in
  // pixels). We can't use the actual number of pixels in radians,
  // as only six pixels  would cause a full 360 degree rotation.
  // So we use a mousesensitivity variable that can be changed to
  // vary how many radians we want to turn in the x-direction for
  // a given mouse movement distance

  // We have to remember that positive rotation is counter-clockwise. 
  // Moving the mouse down is a negative rotation about the x axis
  // Moving the mouse right is a negative rotation about the y axis
  MouseDirection.x = (MiddleX - MouseX)/MouseSensitivity; 
  MouseDirection.y = (MiddleY - MouseY)/MouseSensitivity;

  CurrentRotationX += MouseDirection.y;
  
  // We don't want to rotate up more than one radian, so we cap it.
  if(CurrentRotationX > 1)
  {
    CurrentRotationX = 1;
    return;
  }
  // We don't want to rotate down more than one radian, so we cap it.
  if(CurrentRotationX < -1)
  {
    CurrentRotationX = -1;
    return;
  }
  else
  {
    // get the axis to rotate around the x-axis. 
    Vector Axis = CrossProduct(View - Position, Up);
    // To be able to use the quaternion conjugate, the axis to
    // rotate around must be normalized.
    Axis = Normalize(Axis);

    // Rotate around the y axis
    RotateCamera(MouseDirection.y, Axis.x, Axis.y, Axis.z);
    // Rotate around the x axis
    RotateCamera(MouseDirection.x, 0, 1, 0);
  }
}

void RotateCamera(double Angle, double x, double y, double z)
{
  quaternion temp, quat_view, result;

  temp.x = x * sin(Angle/2);
  temp.y = y * sin(Angle/2);
  temp.z = z * sin(Angle/2);
  temp.w = cos(Angle/2);

  quat_view.x = View.x;
  quat_view.y = View.y;
  quat_view.z = View.z;
  quat_view.w = 0;

  result = mult(mult(temp, quat_view), conjugate(temp));

  View.x = result.x;
  View.y = result.y;
  View.z = result.z;
}

gluLookAt(Position.x, Position.y, Position.z,
          View.x, View.y, View.z, Up.x, Up.y, Up.z).

This document is to be distributed for free and without any modification from its original state. The author declines all responsibility in the damage this document or any of the things you will do with it might do to anyone or to anything. This document and any of its contents is not copyrighted and is free of all rights, you may thus use it, modify it or destroy it without breaking any international law. However according to the author's will, you may not use this document for commercial profit directly, but you may use indirectly its intellectual contents; in which case I would be pleased to receive a mail of notice or even thanks. This is my first tutorial and I am still a student, you must assume that this document is probably not free of small errors and bugs. In the same state of mind, those algorithms are not fully optimised, they are explained for pedagogical purposes and you may find some redundant computations or other voluntary clumsiness. Please be indulgent and self criticise everything you might read. Hopefully, lots of this stuff was taken in sources and books of reference; as for the stuff I did: it has proven some true efficiency in test programs I made and which work as wanted. As said in the introduction: If you have any question or any comment about this text, please send it to the above email address, I'll be happy to answer as soon as possible.

Introduction

Simulating a physical phenomena which obeys to known mathematical equations is, with a number of approximations, always feasable. But what about more abstract concepts, such as feelings, which do not follow any laws? The simplest things we can feel are often the hardest things to capture in a program. Beat detection follows this rule : feeling the beat of a song comes naturally to humans or animals. Indeed it is only a feeling one gets when listening to a melody, a feeling which will make you dance in rhythm or hit a table with your hands on the melody beats. Therefore, how can we teach this beat detection to a machine that can only compute logical operations? In fact there are a number of algorithms which manage to approximate, more or less accurately, this beat detection. We will first study the statistical approach of beat detection on a streaming source and secondly a filtering approach of rhythm extraction on a static song.

This guide assumes the reader has basic signal processing understanding (FFT, convolutions and correlations should sound common) maybe some stuff in statistics will also help (Variance, Average, Principal Components Analysis, will be quoted among others). The point here is not to actually write the code of these algorithms, but more to understand how they work and to be able to adapt or create the appropriate algorithm to a situation. If you have a question or a comment about this text, please send it to the above email address, I'll be happy to answer as soon as possible. Anyway, the aim here is to give more precise ideas on the subject of beat detection to the reader. Enjoy.

I – Statistical streaming beat detection

1 – Simple sound energy

a - A first analysis

The human listening system determines the rhythm of music by detecting a pseudo – periodical succession of beats. The signal which is intercepted by the ear contains a certain energy, this energy is converted into an electrical signal which the brain interprets. Obviously, The more energy the sound transports, the louder the sound will seem. But a sound will be heard as a beat only if his energy is largely superior to the sound's energy history, that is to say if the brain detects a brutal variation in sound energy. Therefore if the ear intercepts a monotonous sound with sometimes big energy peaks it will detect beats, however, if you play a continuous loud sound you will not perceive any beats. Thus, the beats are big variations of sound energy. This first analysis will bring us to our simplest model : Sound energy peaks.

In this model we will detect sound energy variations by computing the average sound energy of the signal and comparing it to the instant sound energy. Lets say we are working in stereo mode with two lists of values : (an) and (bn). (an) contains the list of sound amplitude values captured every Te seconds for the left channel, (bn) the list of sound amplitude values captured every Te seconds for the right channel. So we want to compute the instant energy and the average energy of the signal. The instant energy will in fact be the energy contained in 1024 samples (1024 values of a[n] and b[n]), 1024 samples represent about 5 hundreds of second which is pretty much 'instant'. The average energy should not be computed on the entire song, some songs have both intense passages and more calm parts. The instant energy must be compared to the nearby average energy, for example if a song has an intense ending, the energy contained in this ending shouldn't influence the beat detection at the beginning. We detect a beat only when the energy is superior to a local energy average. Thus we will compute the average energy on say : 44032 samples which is about 1 second, that is to say we will assume that the hearing system only remembers of 1 second of song to detect beat. This 1 second time (44032 samples) is what we could call the human ear energy persistence model; it is a compromise between being to big and taking into account too far away energies, and being too small and becoming to close to the instant energy to make a valuable comparison.

The history buffer where we will keep the last 44032 samples wil contain in fact too lists of samples (B[0]) and (B[1]) corresponding to the left (an) and to the right (bn) channels history.

Simple sound energy algorithm #1:

Every 1024 samples:

• Use the 1024 new samples taken in a[n] and b[n] to compute the instant energy 'e', using the following formula (i0 is the position of the 1024 samples to process):  
  (R1)
• Compute the local average energy '<E>' on the 44100 samples of a history buffer (B):
  (R2)
• Shift the 44032 history buffer (B) of 1024 indexes to the right so that we make room for the 1024 new samples and evacuate the oldest 1024 samples.
• Move the 1024 new samples on top of the history buffer.
• Compare 'e' to 'C * ' where C is a constant which will determine the sensibility of the algorithm to beats. A good value for this constant is 1.3. If 'e' is superior to 'C * ' then we have a beat!

b - Some direct optimisations

This was the basic version of the algorithm, its speed and accurecy can be improved quite easily. The algorithm can be optimised by keeping the energy values computed on 1024 samples in history instead of the samples themselves, so that we don't have to compute the average energy on the 44100 samples buffer (B) but on the instant energies history we will call (E). This sound energy history buffer (E) must correspond to approximately 1 second of music, that is to say it must contain the energy history of 44032 samples (calculated on groups of 1024) if the sample rate is 44100 samples per second. Thus E[0] will contain the newest energy computed on the newest 1024 samples, and E[42] will contain the oldest energy computed on the oldest 1024 samples. We have 43 energy values in history, each computed on 1024 samples which makes 44032 samples energy history, which is equivalent to 1 second in real time. The count is good. The value of 1 second represents the persistance of the music energy in the human ear, it was obtain with experimentations but it may varry a little from a person to another, just adjust it as you feal. So here is what the algorithm becomes:

Simple sound energy algorithm #2:

Every 1024 samples:

• Compute the instant sound energy 'e' on the 1024 new sample values taken in (an) and (bn) using the formula (R1)
• Compute the average local energy <E> with (E) sound energy history buffer:
  (R3)
• Shift the sound energy history buffer (E) of 1 index to the right. We make room for the new energy value and flush the oldest.
• Pile in the new energy value 'e' at E[0].
• Compare 'e' to 'C*<E>'.

c - Sensitivity detection

The imediate draw back of this algorithm is the choice of the 'C' constant. For example in techno and rap music beats are quite intense and precise so 'C' should be quite high (about 1.4); whereas for rock and roll, or hard rock which contains a lot of noise, the beats are more confused and 'C' should be low (about 1.1 or 1.0). There is a way, to make the algorithm determine automatically the good choice for the 'C' constant. We must compute the variance of the energies contained in the energy history buffer (E). This variance, which is nothing but the average of ( Energy Values – Energy average = (E) - <E>), will quantify how marked the beats of the song are and thus will give us a way to compute the value of the 'C' constant. The formula to calculate the variance of the 43 E[i]values is described below (R4). Finally, the greater the variance is the more sensitive the algorithm should be and thus the smaller 'C' will become. We can choose a linear decrease of 'C' with 'V' (the variance) and for example when V  200, C -> 1.0 and when V -> 25, C -> 1.45 (R5). This is our new version of the sound energy beat detection algorithm:

Simple sound energy algorithm #3:

   Every 1024 samples:

• Compute the instant sound energy 'e' on the 1024 new samples taken in (an) and (bn) using the following formula (R1).
• Compute the average local energy <E> with (E) sound energy history buffer using formula (R3).
• Compute the variance 'V' of the energies in (E) using the following formula:
  (R4)
  Compute the 'C' constant using a linear degression of 'C' with 'V', using a linear regression with values (R5):
  (R6)
• Shift the sound energy history buffer (E) of 1 index to the right. We make room for the new energy value and flush the oldest.
• Pile in the new energy value 'e' at E[0].
• Compare 'e' to 'C*[lessthanE>', if superior we have a beat!

   Those three algorithms were tested with several types of music, among others : pop, rock, metal, techno, rap, classical, punk. The fact is the results are quite unpredictable. I will only talk about Simple beat detection algorithm #3 as #2 and #1 are only pedagogical intermediates to get to the #3.

Clearly, the beat detection is very accurate and sounds right with techno and rap, the beats are very precise and the music contains very few noise. The algorithm is quite satisfying for that kind of music and if you aim to use beat detection on techno you can stop reading here, the rest won't change anything to your beat detection. However, even if the improvement of the dynamic 'C' calculations ameliorates things alot, the beat detection on punk, rock and hard rock, is sometimes quite approximate. We can feel it doesn't really get the rythm of the song. Indeed the algorithm detects energy peaks. Sometimes you can hear a drum beat which is sank among other noises and which goes trough the algorithm without being detected as a beat.

To explain this phenomena lets say a guitare and flute make alternatively an amplitude constant note. Each time the first finishes the other starts. The note made by the guitare and the note made by the flute have the same energy but the ear detects a certain rhythm because the notes of the instruments are at different pitch. For our algorithm (who is one might say colorblind) it is just an amplitude constant noise with no energy peaks. This partly explains why the algorithm doesn't detect precisely beats in songs with a lot of instruments playing at different rythms and simultaneously. Our next analysis will make us walk through this difficulty.

Comparing the results we have obtained with the Simple beat detection algorithm #3 to its computing cost, this algorithm is very efficient. If you are not looking for a perfect beat detection than I recommend you use it. Here is a screenshot of a program I made using this algorithm. You fill find the binaries and the sources on my homepage.

Figure 1: The spectrum analyser is not useful for the beat detection it is only for visual matters, but you can see at the top some of the parameters the program computes to execute the algorithm.

2 – Frequency selected sound energy

a - The idea and the algorithm

The issue with our last analysis of beat detection is that it is colorblind. We have seen that this could raise quite a few problems for noisy like songs in rock or pop music. What we must do is give our algorithm the abbility to determine on which frequency subband we have a beat and if it is powerful enough to take it into account. Basically we will try to detect big sound energy variations in particular frequency subbands, just like in our last analysis; unless this time we will be able to seperate beats regarding their color ( frequency subband ). Thus If we want to give more importants to low frequency beats or to high frequency beats it should be more easy. Notice that the energy computed in the time domain is the same as the energy computed in the frequency domain, so we don't have any difference between computing the energy in time domain or in frequency domain. For maths freaks this is called the Parseval Theorem.

Okay that was just a bit of sport, lets go back to the mainstream; Here is how the Frequency selected sound energy algorithm works: The source signals are still coming from (an) and (bn). (an) and (bn) can be taken from a wave file, or directly from a streaming microphone or line input. Each time we have accumulated 1024 new samples, we will pass to the frequency domain with a Fast Fourier Transform (FFT). We will thus obtain a 1024 frequency spectrum. We then divide this spectrum into however many subbands we like, here I will take 32. The more subbands you have, the more sensitive the algorithm will be but the harder it will become to adapt it to lots of different kinds of songs. Then we compute the sound energy contained in each of the subbands and we compare it to the recent energy average corresponding to this subband. If one or more subbands have an energy superior to their average we have detected a beat.

The great progress with the last algorithm is that we now know more about our beats, and thus we can use this information to change an animation, for example. So here is more precisely the Frequency selected sound energy algorithm #1:

   Every 1024 samples:

• Compute the FFT on the 1024 new samples taken in (an) and (bn). The FFT inputs a complex numeric signal. We will say (an) is the real part of the signal and (bn) the imagenary part. Thus the FFT will be made on the 1024 complex values of:
  You can find FFT tutorials and codes in C, Visual Basic or C++ on my homepage in the 'tutorials' section or by typing 'FFT' on Google.
• From the FFT we obtain 1024 complex numbers. We compute the square of their module and store it into a new 1024 buffer. This buffer (B) contains the 1024 frequency amplitudes of our signal.
• Divide the buffer into 32 subbands, compute the energy on each of these subbands and store it at (Es). Thus (Es) will be 32 sized and Es[i] will contain the energy of subband 'i':
  (R7)
• Now, to each subband 'i' corresponds an energy history buffer called (Ei). This buffer contains the last 43 energy computations for the 'i' subband. We compute the average energy <E> for the 'i' subband simply by using:
  (R8)
• Shift the sound energy history buffers (Ei) of 1 index to the right. We make room for the new energy value of the subband 'i' and flush the oldest.
• Pile in the new energy value of subband 'i' : Es[i] at Ei[0].
  (R9)
• For each subband 'i' if Es[i] > (C*<E>) we have a beat !

   To help out visualising how the data piles work have a look at this scheme:

Figure 2: This is how the energy data is organized.

Now the 'C' constant of this algorithm has nothing to do with the 'C' of the first algorithm, because we deal here with separated subbands the energy varies globally much more than with colorblind algorithms. Thus 'C' must be about 250. The results of this algorithm are convincing, it detects for example a symbal rhythm among other heavy noises in metal rock, and indeed the algorithm separates the signal into subbands, therefore the symbal rhythm cannot pass trough the algorithm without being recognized because it is isolated in the frequency domain from other sounds. However the complexity of the algorithm makes it useful only if you are dealing with very noisy sounds, in other cases, Simple beat detection algorithm #3 will do the job.

b - Enhancements and beat decision factors.

There are ways to enhance a bit more our Frequency selected sound energy algorithm #1.

First we will increase the number of subbands from 32 to 64. This will take obviously more computing time but it will also give us more precision in our beat detection. The second way to develop the accuracy of the algorithm uses the defaults of human ears. Human hearing system is not perfect; in fact its transfer function is more like a low pass filter. We hear more easily and more clearly low pitched noises than high pitch noises. This is why it is preferable to make a logarithmic repartition of the subbands. That is to say that subband 0 will contain only say 2 frequencies whereas the last subband, will contain say 20. More precisely the width 'wi' of the 'n' subbands indexed 'i' can be obtained using this argument:

• Linear increase of the width of the subband with its index:
  (R10)
• We can choose for example the width of the first subband:
  (R11)
• The sum of all the widths must not exceed 1024 ((B)'s size):
  (R12)

   Once you have equations (R11) and (R12) it is fairly easy to extract 'a' and 'b', and thus to find the law of the 'wi'. This calculus of 'a' and 'b' must be made manually and 'a' and 'b' defined as constants in the source; indeed they do not vary during the song.

So in fact in Frequency selected sound energy algorithm #1, all we have to modify is the number of subbands we will take equal to 64 and the (R7) relation. This relation becomes:

(R7)'

It may seem rather complicated but in fact it is not. Replacing this relation (R7) with (R7)' we have created Frequency selected sound energy algorithm #2. If you have musics with very tight and rapid beats, you may want to compute the stuff more frequently than every 1024 samples, but this is only for special cases, normally the beat should not be shorter than 1/40 of second.

Using the advantages of Frequency selected beat detection you can also enhance the beat decision factor. Up to now it was based on a simple comparison between the instant energy of the subband and the average energy of the subband. This algorithm enables you to decide beats differently. You may want for examples to cut beats which correspond to high pitch sounds if you run techno music or you may want to keep only [50-4000Hz] beats if you are working with speech signal. This algorithm has the advantage of being perfectly adaptable to any kind or category of signal which was not the case of Simple beat detection algorithm #3. Notice that the correspondants between index 'i' of the FFT transform and real frequency is given by formula:

• If 'i' < 'N/2' then:
  (R13)
• Else 'i' >= 'N/2' then:
  (R14)

   So 'i' is the index of the value in the FFT output buffer, N is the size of the FFT transform (here 1024), fe is the sample frequency (here 44100). Thus index 256 corresponds to 10025 Hz. This formula may be useful if you want to create your own subbands and you want to know what the correspondants between indexes and real frequency are.

Another way of filtering the beats, or selecting them, is choosing only those which are marked and precise enough. As we have seen before, to detect the accuracy of beats we must compute the variance of the energy histories for each subband. If this variance is high it means that we have great differences of energy and thus that the beat are very intense. Thus all we have to do is compute this variance for each subband and add a test in the beat detection decision. To the "Es[i] > (C*<Ei>)" condition we will add "and V((Ei))>V0". V0 will be the variance limit, with experience 150 is a reasonable value. Now the V((Ei)) value is easy to compute, just follow the following equality if you don't see how:

(R15)

The last (finally) way to enhance your beat detection, is to make the source signal pass through a derivation filter. Indeed differentiating the signal makes big variations of amplitude more marked and thus makes energy variations more marked and recognisable later in the algorithm. I haven't tried this optimisation but according to some sources this is quite useful. If you try it please give me your opinion on it!

Concerning the results of the Frequency selected sound energy algorithm #2 I must admit they are way more satisfying than the Simple sound energy algorithm #3. In a song the algorithm catches the bass beats as well as the tight cymbals hits. I insist on the fact that you may also select the beats in very different ways, which becomes quite useful if you know you are going to run techno music with low pitch beats for example. You may also select the beats differently according to there accuracy with the variance criteria. There are many other ways to decide beats; it is up to you to explore them and find the one which fits the most your needs.

I used Frequency selected sound energy algorithm #2 algorithm in a demo program of which you can see some screenshots just below. One can see quite clearly that there is a beat in low frequencies (probably a bass or drum hit) and also a high pitch beat (probably a cymbal or such):

Figure 3: The histogram represents the instant energy of the 128 subbands.

Figure 4: On this screenshot you can see the 128 subbands E/<E> ratios, and the bandwidth of the subbands. When the algorithm detects a beat in a subband, it overwrites the ratio in red.

The reason why I called this part of the document, 'Statistical streaming beat detection', is that the energy can be considered as a random variable, of which we have been calculating values over time. Thus we could interpret those values as a sampling for the statistical analysis of a random variable. But one can push this approach further. When we have separated the energies histories into 128 subbands we created in fact 128 random energy variables. We can then apply some of the general statistical methods of analysis. For example, the principal components analysis method will enable you to determine if some of the subbands are directly linked or independent. This would help us to regroup some subbands which are directly linked and thus make the beat decision more efficient. However, this method is basically just far too computing expensive and maybe just too hard to implement comparing to the results we want. If you are looking for a really good challenge in beat detection you could push in this direction.

II – Filtering rhythm detection

1 – Derivation and Comb filters

a - Intercorrelation and train of impulses

While in the first part of this document we had seen beat detection as a statistical sound energy problem, we will approach it here as a signal processing issue. I haven't tried myself this algorithm but I greatly inspired myself of some work which was done and tested with this algorithm (see the sources section), so it really works but I won't detail so much its implementation as I did in the first part. The basic idea is the following. If we have two signals x(t) and y(t), we can evaluate a function which will give us an idea of how much those two signals are similar. This function called 'inter correlation function' is the following:

(R16)

For our purposes we will always take alpha=1. This function quantifies the energy that can be exchanged between the signal x(t) and the signal y(t – b), it gives an evaluation of how much those two functions look like each other. The role of the 'b' is to eliminate the time origin problem; two signals should be compared without regarding their possible time axis translation. As you may have noticed, we are creating this algorithm for a computer to execute; basically we have to re-write this formula in discrete mode, it becomes:

(R17)

Now the point is we can valuate the similarity of our song and of a train of impulses using this function. The train of impulses being at a known frequency (the beats per minute we want to test) is what we call a comb filter. The energy of the Y function gives us an evaluation of the similarity between our song and this train of impulses, thus it quantifies how much the rhythm of the train of impulses is present in the song. If we compute the energy of the ? function for different train of impulses, we will be able to see for which of them the energy is the biggest, and thus we will know what is the principal Betas Per Minute rate. This is called combfilter processing. We will note 'Ey' this energy, where x[k] is our song signal and y[k] our train of impulses:

(R18)

By the way, here is what a train of impulses looks like:

Figure 5: This is the train of impulses function, it is caracterised by its period Ti.

So the period Ti of the impulses must correspond to the beats per minutes we want to test on our song. The formula that links a beat per minute value with the Ti in discrete time is the following:

(R19)

fs is the sample frequency, if you are using good quality wave files, it is usually of 44100. BPM is the beats per minute rate. Finally, Ti is the number of indexes between each impulse.

   Now because it is quite computing expensive to compute the (R17) formula, we pass to the frequency domain with a FFT and compute the product of the X[v] and Y[v] (FFT's of x[k] and y[k]). Then we compute the energy of the ? function directly in the frequency domain, thus Ey is given with by the following formula:

(R17)'

b - The algorithm

Note that this algorithm is much too computing expensive to be ran on a streaming source, or on a song in a whole. It is executed on a part of the song, like 5 seconds taken somewhere in the middle. Thus we assume that the tempo of the song is overall constant and that the middle of the song is tempo-representative of the rest of the song. So finally here are the steps of our algorithm:

Derivation and Combfilter algorithm #1:

• Choose roughly 5 seconds of data in the middle of the song and copy it to the 'a[k]' and 'b[k]' lists. The dimension of those lists is noted N. As usual 'a[k]' is the array for the left values, and 'b[k]' is the array for the right values.
• Compute the FFT of the complex signal made with a[k] for the real part of the complex signal and b[k] the imaginary part of the complex signal as seen in Frequency selected sound energy algorithm #1. Store the real part of the result in ta[k] and the imaginary in tb[k].
• For all the beats per minute you want to test, for example from 60 BPM to 180 BPM per step of 10, we will note BPMc the current BPM tested:
  
  • Compute the Ti value corresponding to BPMc using formula (R19).
  • Compute the train of impulses signal and store it in (l) and (j). (l) and (j) are purely identical, we take two lists of values to have a stereo signal so that we don't have dimension issues further on. Here is how the train of impulses is generated:
    
    for (k = 0; k < N; k++)
    {
      if ((k % Ti) == 0)
        l[k] = j[k] = AmpMax;
      else
        l[k] = j[k] = 0;
    }
  • Compute the FFT of the complex signal made of (l) and (j). Store the result in tl[k] and tj[k].
  • Finally compute the energy of the correlation between the train of impulses and the 5 seconds signal using (R17)' which becomes:
    (R20)
  • Save E(BPMc) in a list.
  
  


• The rhythm of the song is given by the BPMmax, where the max of all the E(BPMc) is taken. We have the beat rate of the song!

   The AmpMax constant which appears in the train of impulse generation, is given by the sample size of the source file. Usually the sample size is 16 bits, for a good quality .wav. If you are working in non-signed mode, the AmpMax value will be of 65535, and more often if you work with 16 bits signed samples the AmpMax will be 32767.

One of the first ameliorations that could be given to Derivation and Combfilter algorithm #1 is indeed adding a derivation filter before the combfilter processing. To make beats more detectable, we differentiate the signal. This accentuates when the sound amplitude changes. So instead of dealing with a[k] and b[k] directly we first transform them using the formula below. This modification added before the second step of Derivation and Combfilter algorithm #1 constitutes Derivation and Combfilter algorithm #2.

(R21)

As I said before, I haven't tested this algorithm myself; I can't compare its results with the algorithms of the first part. However, throughout the researches I made this seems to be the algorithm universally accepted as being the most accurate. It does have disadvantages: Great computing power consumption, can only be done on a small part of a song, assumes the rhythm is constant throughout the song, no possible streaming input (for example from a microphone), the signal needs to be known in advance. However it returns a deeper analysis than the first part algorithms, indeed for each BPM rate you have a value of its importance in the song.

2 – Frequency selected processing

a - Curing colorblindness

As in the first part of the tutorial, we have the same problem our last algorithm doesn't make the difference between cymbals beat or a drum beat, so the beat detection becomes a bit dodgy on very noisy music like hard rock. To heal our algorithm we will proceed as we had done in Frequency selected sound energy algorithm #1 we will separate our source signal into several subbands. Only this time we have much more computing to do on each subband so we will be much more limited in their number. I think that 16 is a good value. We will use a logarithmic repartition of these subbands as we had done before.

Let's modify the Derivation and Combfilter algorithm #2. The values particular to a subband will be characterized with a little 's' at the end of the name of the variable. So we will separate ta[k] and tb[k] into 16 subbands. Each subbands values array will be called tas[k] and tbs[k] ('s' varies from 1 to 16).

Frequency selected processing combfilters algorithm #1:

• Choose roughly 5 seconds of data in the middle of the song and copy it to the 'a[k]' and 'b[k]' lists. The dimension of those lists is noted N. As usual 'a[k]' is the array for the left values, and 'b[k]' is the array for the right values.
• Differentiate the a[k] and b[k] signals using (R21).
• Compute the FFT of the complex signal made with a[k] for the real part of the complex signal and b[k] the imaginary part of the complex signal as seen in Frequency selected sound energy algorithm #1. Store the real part of the result in ta[k] and the imaginary in tb[k].
• Generate the 16 subband array values tas[k] and tbs[k] by cutting ta[k] and tb[k] with a logarithmic rule.
• For all the subbands tas[k] and tbs[k] (s goes from 1 to 16). Ws is the length of the 's' subband.
• For all the beats per minute you want to test, for example from 60 BPM to 180 BPM per step of 10, we will note BPMc the current BPM tested :
  
  • Compute the Ti value corresponding to BPMc using formula (R19).
  • Compute the train of impulses signal and store it in (l) and (j). (l) and (j) are purely identical, we take two lists of values to have a stereo signal so that we don't have dimension issues further on. Here is how the train of impulses is generated:
    
    for (k = 0; k < ws; k++)
    {
      if ((k % Ti) == 0)  
        l[k] = k[k] = AmpMax;
      else  
        l[k] = j[k] = 0;
    }
  • Compute the FFT of the complex signal made of (l) and (j). Store the result in tl[k] and tj[k].
    Finally compute the energy of the correlation between the train of impulses and the 5 seconds signal using (R17)' which becomes:
    (R22)
  • Save E(BPMc,s) in a list.
  
  


• The rhythm of the subband is given by the BPMmaxs, where the max of all the E(BPMc,s) is taken (s is fixed). We will call this max EBPMmaxs. We have the beat rate of the subband we store it in a list.


• We do this for all the subbands, we have a list of BPMmaxs for all subbands, we can than decide which one to consider.

   As in Frequency selected sound energy algorithm #2 we will then be able to decide the rhythm according to frequency bandwidth criteria. Or, if you want to take all the subbands into account you can compute the final BPM of the song, by calculating the barycentre of all the BPMmaxs affected with their max of E(BPMc,s). Like this:

(R23)

I must admit I haven't concretely tested this algorithm. Others have done this already and here an overview of the results for Derivation and Combfilter algorithm #2 (source: http://www.owlnet.ri...c/beatalgo.html).

Figure 6: This is plot of the EBPMc values function of BPMc. The algorithm will take the max of the EBPMc value to find the final BPM of the song. Here we see this max is reached for BPMc=75. The final BPM is 75!

Conclusion

Finding algorithms for beat detection is very frustrating. It seems so obvious to us, humans, to hear those beats and somehow so hard to formalise it. We managed to approximate more or less accurately and more or less efficiently this beat detection. But the best algorithms are always the ones you make yourself, the ones which are adapted to your problem. The more the situation is precise and defined the easier it is! This guide should be used as a source of ideas. Even if beat detection is far from being a crucial topic in the programming scene, it has the advantage of using lots of signal processing and mathematical concepts. More than an end in itself, for me, beat detection was a way to train to signal processing. I hope you will find some of this stuff useful.

Sources and Links

• http://www.owlnet.ri...c/beatalgo.html: This site greatly inspired me for the second part of the tutorial, very well explained.
• http://www.cs.prince...rs/amta2001.pdf: Audio analysis using the discrete wavelet transform.
• Tempo and beat analysis of acoustic musical signals by Eric D. Scheirer
• Pacemaker by Olli Parviainen

• Translation
• Scale
• Rotation