How To Load/Save Data in Python

Loading/Saving data is a common task in numerical analysis. In python, the most convenient way is to use the save/load function provided in numpy package. An example is as the follows

import numpy as np
rand_array = np.random.rand(1,10);
np.save('rand.npy',rand_array);
b=np.load('rand.npy');

This example shows how to save an array in a .npy file and then load. This save/load function is quite similar to the save/load in Matlab except for that Matlab allows saving multiple arrays in single file while in Python one file can hold only one array. For comparison purpose, Matlab save/load code is shown below where two arrays with random numbers are save in the same file

rand_array=rand(1,10);
rand_array2 = rand(1,10);
save('rand_file','rand_array','rand_array2');
myFile = load('rand_file');

3D LIDAR DIY Using PulsedLight LIDAR-Lite Device and Servo

LIDAR stands for LIght Detection And Ranging. LIDAR uses laser to measure the distance. LIDAR is capable to do things no many other sensors can do. It can accurately measure distance (tens to hundreds of meters) of objects around, which generates a good picture of surrounding environment. Because of this functionality, LIDAR is a must-have for many self-driving car designs. Many tech hobbyists want to put their fingers on a LIDAR device. However, commercial LIDAR devices are quite expensive and they are usually well above one thousand dollars. Good news is that with availability of cheap laser component like PulsedLight LIDAR-Lite and open source platform like Arduino, now we can make by ourselves a 3D LIDAR scanner with cost well below 1000$. In this blog, I will introduce such a device made by myself.

PulsedLight LIDAR-Lite is used in this device to measure the distance using laser. The maximum distance LIDAR-Lite measures is 40 meters. Two servos, one for vertical direction and another for horizontal direction, are used for scanning the 3D space. Both horizontal and vertical servos can rotate up to 180 degrees. This defines the coverage of this scanner. Arduino platform controls all of them. By connecting Arduino to PC, 3D scan data can be downloaded to PC to control point cloud file. At the same time, PC can power the device through USB. Thus the device does not need to have its own power. Figure below shows how the device looks like.

























When the device scans, it first fixes vertical angle and sweeps horizontally. Then it moves one vertical step and sweeps horizontally again. This Youtube video shows the scanning process

After scanning finished, a point cloud will be generated. Often several point clouds need to be stitched together to generate the whole picture for a building. The residential house model below is generated in that way.

3D model can be produced based on point cloud (Thanks Igor!). Below is the 3D model for the same residential house.

After 3D printing, we have a plastic model of the house. We can view it together with Google street map side-by-side.

Signal Processing Magic (3) -- Farrow Structure

Farrow structure was invented by C. W Farrow. It minimizes the number of multiplication used in fractional interpolation. Farrow structure is often deployed when there is fractional sampling rate change in the signal path. For example, ADC rate is set as 120MHz but baseband signal processing only supports 100MHz. Therefore, a filter needs to be designed for supporting rate conversion from 120MHz to 100MHz. Farrow structure is one candidate for this kind of job.

The setting of Farrow depends on the polynomial it supports. Let us take an example of using cubic Lagrange for fractional interpolation. The output y[n] can be written as combination of input x[n] as:
\[y[n]=C_{-2}x[n-2]+C_{-1}x[n-1]+C_{0}x[n]+C_{1}x[n+1]\]

Assuming fractional delay is \(\tau\), using Lagrange formula, \(C_{-2}\) equals to the division of \((-\tau-(-1))(-\tau)(-\tau-1)\) by \((-2-(-1))(-2)(-2-1)\). Thus,
\[C_{-2}=\frac{\tau^{3}}{6}-\frac{\tau}{6}\]
Extend the same formula to other coefficients, we have
\[C_{-1}=-\frac{\tau^{3}}{2}+\frac{\tau^{2}}{2}+\tau\]
\[C_{0}=\frac{\tau^{3}}{2}-\tau^{2}-\frac{\tau}{2}+1\]
\[C_{1}=-\frac{\tau^{3}}{6}+\frac{\tau^{2}}{2}-\frac{\tau}{3}\]
Rewriting the y[n] equation, we have
\[y[n]=\frac{1}{6}[(x[n-2]-3x[n-1]+3x[n]-x[n+1])\tau^{3}+ \\
(3x[n-1]-6x[n]+3x[n+1])\tau^2+ \\
(-x[n-2]+6x[n-1]-3x[n]-2x[n+1])\tau+6x[n]] \]

Corresponding Farrow structure is shown below. It only needs three multipliers and \(\tau\) is input to all multipliers. This simplifies the hardware design. When polynomial is different, the coefficients of FIR filter change accordingly. We should also note that due to the difference between input and output rates, one input sample not necessarily generates one output sample. A control logic is needed to dictate when to generate outputs.

Signal Processing Magic (2) -- Sigma Delta ADC

Sigma Delta ADC is a commonly used architecture for converting analog signal to digital. The conception of ADC is quite simple: let us say we have 4 bits with bit3 representing 1volt, bit2 for 0.5 v, bit1 for 0.25v and bit0 for 0.125v, I can quantize an analog input in volt to these 4 bits. For example, if the input is 1.2 volts, then it becomes [bit3 bit2 bit1 bit0] = [1010] = 1.25v. The quantization error is equal to 1.2-1.25 = -0.05v. Quantization noise comes from quantization error. A straightforward thought is to use binary search to convert analog input to digital. This is the basic idea behind SAR (Successive Approximation Register) ADC. 

Sigma Delta ADC, our main topic in this article, is constructed differently. Assuming x(n) is the input and y(n) is the output of sigma delta modulator, y(n) has a much higher sampling rate than x(n) and y(n) is also only 1 bit. So in term of y(n), sigma delta can be seen as a trade off between sampling rate and bitwidth. y(n) will then pass through lowpass filter and decimator to produce the multiple-bit ADC output we commonly work with.

Now let us do some math to figure out what really happens in sigma delta. First we calculate the transfer function between v(n) and u(n)

\[z^{-1}(U(z)+V(z))=V(z)\]
which gives
\[H(z)=\frac{V(z)}{U(z)}=\frac{z^{-1}}{1-z^{-1}}\]

The transfer function from input, x(n), to output, y(n), is equal to
\[\frac{Y(z)}{X(z)}=\frac{H(z)}{1+H(z)}=z^{-1}\]
which tell us there is signal delay but no signal distortion.

e(n) added to 1-bit ADC in the diagram represents quantization noise. The transfer function from e(n) to y(n) is
\[\frac{Y(z)}{E(z)}=\frac{1}{1+H(z)}=1-z^{-1}\]
Frequency response of Y(z)/E(z) is shown below. Its amplitude increases with frequency which means sigma delta modulator pushes out noise from in-band to out-of-band. Reducing in-band noise level is the major benefit of sigma delta. What we discussed here is first-order sigma delta. If we goes to higher order sigma delta, the in-band noise level will be further reduced but the cost paid is more sophisticated hardware and higher out-of-band noise level.

Signal Processing Magic (1) -- CIC Filter

In the world of signal processing, there are many elegant schemes. They are beautiful and pragmatic at the same time. Sometime you are wondering how the inventors create them. I feel obliged to introduce them to my audience. Today, let me start from a filter named CIC (Cascaded Integrator–Comb). CIC filter is mainly used for integer up-sampling and down-sampling. The main benefit of this scheme is that no multiplication is needed, which is welcomed in term of hardware complexity.

Figure below shows CIC filter architecture with N stages and decimation rate of R. M is an extra parameter and it can be 1 or 2. Normally M is set as 1. CIC filter can also be used as interpolation and the flow is opposite when used as decimation.

Accordingly, we can write the transfer function of decimation CIC filter as:

\[(\frac{1-z^{-RM}}{1-z^{-1}})^{N}\]


By increasing N with fixed R, the filter can have more rejection.

By increasing R with fixed N, the filter's cutoff frequency decreases proportionally.

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.

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!