Saturday, October 21, 2017

Digest of ESPRIT Algorithm

For most folks, ESPRIT means a fashion branch (http://www.esprit.eu/). But for nerds, ESPRIT stands for Estimation of Signal Parameter via Rotational Invariance Technique. It is a signal processing algorithm used to process things like signal coming from array. ESPRIT is often mentioned together with MUSIC algorithm. Again, MUSIC is not the real music played by Mozart but stands for MUltiple SIgnal Classification.

Recently I spent some time to understand the details of ESPRIT algorithm and here is what I have learnt. The biggest difference between ESPRIT and MUSIC is that ESPRIT has dual receivers. The diagram below is from a tutorial found in Internet (http://www.girdsystems.com/pdf/GIRD_Systems_Intro_to_MUSIC_ESPRIT.pdf). As shown here, to make ESPRIT algorithm possible, there are many pairs of sensors. The angle of arrive signal is obtained by comparing these two received signals.



















The best way for understanding such algorithm is to run a Matlab script. The website of a book named "Spectral Analysis of Signals" provide a well written example of Matlab code for ESPRIT. The code is as below:


function w=esprit(y,n,m)
%
% The ESPRIT method for frequency estimation.
%
%  w=esprit(y,n,m);
%
%      y  ->  the data vector
%      n  ->  the model order
%      m  ->  the order of the covariance matrix in (4.5.14)
%      w  <-  the frequency estimates
%

% Copyright 1996 by R. Moses

y=y(:);
N=length(y);                       % data length

% compute the sample covariance matrix
R=zeros(m,m);
for i = m : N,
   R=R+y(i:-1:i-m+1)*y(i:-1:i-m+1)'/N;
end

% to use the forward-backward approach, uncomment the next line
% R=(R+fliplr(eye(m))*R.'*fliplr(eye(m)))/2;

% get the eigendecomposition of R; use svd because it sorts eigenvalues
[U,D,V]=svd(R);
S=U(:,1:n);

phi = S(1:m-1,:)\S(2:m,:);

w=-angle(eig(phi));
return

R in the code is a m x m matrix. The algorithm divides S(1:m, :) to two parts of S(1:m-1, :) and S(2:m, :). It is like dividing a large array to two smaller arrays with overlap of m-2 elements. I think this can be further fine tuned. For example, S(1:m, :) can be divided to S(1:m-2, :) and S(3:m, :) with one less overlapping element. "S(1:m-1,:)\S(2:m,:)" is to find least square solution of S(1:m-1,:)*X = S(2:m,:). By this method, it finds the delta between received signal of these two slightly different arrays.

Sunday, October 15, 2017

Image Capture Using Arduino UNO and ArduCam OV2640 Module

Today we run through an example of how to connect a camera module to Arduino and get image captured. 

Step 1: I already had a Arduino Uno board. Thus the first step is to find a camera module which can be connected to this board. After some googling, I found OV2640 module made by ArduCam. This module has a 2MP image sensor and can be conveniently connected to Uno board.

Step 2: Connecting HW. Below is a Arduino Uno pinout diagram (from http://www.electroschematics.com/7958/arduino-uno-pinout/). ArduCam OV2640 module has eight pins. Seven out of these eight pins are connected to the upper left corner of Uno board circled in green in the diagram.



As the figure below show, except for the pin of 5v, which is connected to the 5v pinout in lower end of the board, all other seven pins are connected to the upper left corner.


Step 3: SW setup. For code to be loaded to Uno board, we use an example program provided by ArduCam named ArduCAM_SPI_OV2640_FIFO_UART.ino. The link to this program from ArduCam's website is broken as on 2017/10/25. But the program can still be found in https://github.com/sumotoy/ArduCAM/tree/master/examples/REVC/ArduCAM_SPI_OV2640_FIFO_UART

ArduCAM_SPI_OV2640_FIFO_UART.ino enables camera configuration setup and image capture. In order to see the image from PC, we still need a host program. For host program in PC, we can use an ArduCam app named ArduCAM_host_v1.0.exe which can be downloaded from https://github.com/sumotoy/ArduCAM/tree/master/examples/mini

After opening this app from PC, you need to set COM port and Baud rate to match with what is used. Baud rate in default Arduino program is 115200. You need to find out COM port by yourself and Arduino GUI can be a good tool. The image resolution is set as 320x240, which already matches with the default. After setup done, press "Open" button. Then press "Capture" botton. If everything is set up correctly, you should see the image. Bingo!








Sunday, May 7, 2017

Resizable Window for Pygame + OpenGL

Pygame + OpenGL provides a good platform for implementing OpenGL in Python language. It is straightforward to make the OpenGL window resizable. In the code below, we first detect the event of pygame video size change. When this event occurs, the display is reset with the width and height of the new window. The OpenGL perspective is also changed accordingly.

            if event.type == pygame.VIDEORESIZE:
                # For window size change
                surface = pygame.display.set_mode((event.w, event.h), DOUBLEBUF|OPENGL|RESIZABLE)
                gluPerspective(45, (1.0*event.w/event.h), 0.1, 50.0)
                glTranslatef(0.0,0.0, -5)

A sample Python program is listed below which shows a cube in OpenGL window which can be resized:


import pygame
from pygame.locals import *
import math

from OpenGL.GL import *
from OpenGL.GLU import *

verticies = (
    (1, -1, -1),
    (1, 1, -1),
    (-1, 1, -1),
    (-1, -1, -1),
    (1, -1, 1),
    (1, 1, 1),
    (-1, -1, 1),
    (-1, 1, 1)
    )

edges = (
    (0,1),
    (0,3),
    (0,4),
    (2,1),
    (2,3),
    (2,7),
    (6,3),
    (6,4),
    (6,7),
    (5,1),
    (5,4),
    (5,7)
    )


def Cube():
    glBegin(GL_LINES)
    for edge in edges:
        for vertex in edge:
            glVertex3fv(verticies[vertex])
    glEnd()
    #print(repr(verticies[0][0])+','+repr(verticies[0][1]));
 
    glBegin( GL_POINTS );
    glColor3f(1,1,1); 
    for i in range(0,8):
         glVertex3f(verticies[i][0], verticies[i][1], verticies[i][2]);
    glEnd();


def main():
    pygame.init()
 
    display = (1000,750)
    surface = pygame.display.set_mode(display, DOUBLEBUF|OPENGL|RESIZABLE)
    gluPerspective(45, (1.0*display[0]/display[1]), 0.1, 50.0)
    glTranslatef(0.0,0.0, -5)

    while True:
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                pygame.quit()
                quit()
            if event.type == pygame.VIDEORESIZE:
                # For window size change
                surface = pygame.display.set_mode((event.w, event.h), DOUBLEBUF|OPENGL|RESIZABLE)
                gluPerspective(45, (1.0*event.w/event.h), 0.1, 50.0)
                glTranslatef(0.0,0.0, -5)

        glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT)
        Cube()

        pygame.display.flip()
        pygame.time.wait(10)


main()

Thursday, April 13, 2017

Using Mouse for Object Zoom In/Zoom Out/Rotation in Python + OpenGL

OpenGL is supported in multiple languages including Python. A common task performed in OpenGL is to zoom in, zoom out and rotate the 3D object. In this blog, we provide an example of how to carry out these tasks. We utilize a Python package name pyname in our code.

In our code, when one rolls up the wheel, the object will be zoomed in; when one rolls down the wheel, the object will be zoomed out. Zoom in/out is implemented by the function of glScaled. When the scale is larger than 1, the object will be enlarged; when the scale is smaller than 1, the object will be shrink. The code snippet for zoom function is:


    if event.type == pygame.MOUSEBUTTONDOWN and event.button == 4: # wheel rolled up
        glScaled(1.05, 1.05, 1.05);
    elif event.type == pygame.MOUSEBUTTONDOWN and event.button == 5: # wheel rolled down
        glScaled(0.95, 0.95, 0.95);


Rotation function is designed in this way: when you move the mouse  to a direction while pressing the left button, the object rotates to that direction. Since the movement of mouse in the monitor is 2D, the first step is to calculate the mouse movement in the x-axis and y-axis directions. pygame allows to record the 2D location of the mouse. The delta between the current location and the previous location gives dx and dy, which indicates the horizontal and vertical directions of movement:
        x, y = event.pos;
        dx = x - lastPosX;
        dy = y - lastPosY;

As the next step of implementing rotation, model view matrix needs to be retrieved. Since the 3D object keeps rotating, model view matrix tells the current orientation.

            modelView = (GLfloat * 16)()
            mvm = glGetFloatv(GL_MODELVIEW_MATRIX, modelView)

The model view matrix is a 4x4 matrix. The first 3 rows by 3 columns are the rotation matrix which is also a unitary matrix. This rotation matrix, Q, tells the orientation of the object. In order to achieve the desired rotation effect, elements of the rotation matrix need to be multiplied with dx and dy. Assuming Q = [q0|q1|q2], since dy represents the rotation around x-axis and dx represents the rotation around y-axis, q0*dy+q1*dx is combined rotation vector. sqrt(dx^2+dy^2) is the magnitude of rotation. The code is as below:

            temp = (GLfloat * 3)();
            temp[0] = modelView[0]*dy + modelView[1]*dx;
            temp[1] = modelView[4]*dy + modelView[5]*dx;
            temp[2] = modelView[8]*dy + modelView[9]*dx;
            norm_xy = math.sqrt(temp[0]*temp[0] + temp[1]*temp[1] + temp[2]*temp[2]);
            glRotatef(math.sqrt(dx*dx+dy*dy), temp[0]/norm_xy, temp[1]/norm_xy, temp[2]/norm_xy);


For reference, an example Python program is provided below. To run this program, the user needs to install pygame Python package.

import pygame
from pygame.locals import *
import math

from OpenGL.GL import *
from OpenGL.GLU import *

lastPosX = 0;
lastPosY = 0;
zoomScale = 1.0;
dataL = 0;
xRot = 0;
yRot = 0;
zRot = 0;

verticies = (
    (1, -1, -1),
    (1, 1, -1),
    (-1, 1, -1),
    (-1, -1, -1),
    (1, -1, 1),
    (1, 1, 1),
    (-1, -1, 1),
    (-1, 1, 1)
    )

edges = (
    (0,1),
    (0,3),
    (0,4),
    (2,1),
    (2,3),
    (2,7),
    (6,3),
    (6,4),
    (6,7),
    (5,1),
    (5,4),
    (5,7)
    )


def Cube():
    glBegin(GL_LINES)
    for edge in edges:
        for vertex in edge:
            glVertex3fv(verticies[vertex])
    glEnd()
    #print(repr(verticies[0][0])+','+repr(verticies[0][1]));
 
    glBegin( GL_POINTS );
    glColor3f(1,1,1); 
    for i in range(0,8):
         glVertex3f(verticies[i][0], verticies[i][1], verticies[i][2]);
    glEnd();
 
def mouseMove(event):
    global lastPosX, lastPosY, zoomScale, xRot, yRot, zRot;
 
    if event.type == pygame.MOUSEBUTTONDOWN and event.button == 4: # wheel rolled up
        glScaled(1.05, 1.05, 1.05);
    elif event.type == pygame.MOUSEBUTTONDOWN and event.button == 5: # wheel rolled down
        glScaled(0.95, 0.95, 0.95);
 
    if event.type == pygame.MOUSEMOTION:
        x, y = event.pos;
        dx = x - lastPosX;
        dy = y - lastPosY;
        
        mouseState = pygame.mouse.get_pressed();
        if mouseState[0]:

            modelView = (GLfloat * 16)()
            mvm = glGetFloatv(GL_MODELVIEW_MATRIX, modelView)
   
   # To combine x-axis and y-axis rotation
            temp = (GLfloat * 3)();
            temp[0] = modelView[0]*dy + modelView[1]*dx;
            temp[1] = modelView[4]*dy + modelView[5]*dx;
            temp[2] = modelView[8]*dy + modelView[9]*dx;
            norm_xy = math.sqrt(temp[0]*temp[0] + temp[1]*temp[1] + temp[2]*temp[2]);
            glRotatef(math.sqrt(dx*dx+dy*dy), temp[0]/norm_xy, temp[1]/norm_xy, temp[2]/norm_xy);

        lastPosX = x;
        lastPosY = y;
        


def main():
    pygame.init()
 
    display = (1000,750)
    pygame.display.set_mode(display, DOUBLEBUF|OPENGL, RESIZABLE)

    gluPerspective(45, (1.0*display[0]/display[1]), 0.1, 50.0)
    glTranslatef(0.0,0.0, -5)
 

    while True:
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                pygame.quit()
                quit()
            mouseMove(event);

        glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT)
        Cube()
        pygame.display.flip()
        pygame.time.wait(10)


main()

How to Run N Nearest Neighbor Search in Python using kdtree Package

Given S points scattered in a K-dimension space, N nearest neighbor search algorithm finds out for certain point, which N out of S points are its closest neighbors. To implement N nearest neighbor searching algorithm, a kd tree needs to be constructed for all these S points. Searching in the internet, the most popular way to set up kd tree in Python is to use the scipy package. However, I keeps failing on installing scipy package due to the lack of Lapack package. Lapack is a linear algorithm library. Instead, I find another way. Probably some other folks also have the issue of scipy installation. Therefore, let me what I learnt.

My solution is to install the kdtree package. The package can be found in https://github.com/stefankoegl/kdtree. It is pretty straightforward to use. Let me share my example code of N nearest neighbor search algorithm. All points here are in a 3D dimension with N = 2. The average distance to its nearest neighbors is printed out for each point.



import kdtree
import math
import numpy

data = ((1, 0, 0),(2, 0, 0), (2, 0, 0), (3, 0, 0));

# Set up kd tree
myTree = kdtree.create(dimensions=3)
for i in range(0, len(data)):
  tempData = data[i];
  myTree.add((float(tempData[0]), float(tempData[1]), float(tempData[2])));

# Nearest neighbor search
minDist = []
NN = 2; # Number of neighbors to search
for i in range(0, len(data)):
   searchResult = myTree.search_knn(data[i], NN+1); # Find two closest neighbors including itself
   sumDist = 0;
   for p in range(0, NN):
     sumDist = sumDist + math.sqrt(searchResult[p+1][1]);
   avgDist = sumDist/NN;
   print avgDist