Myanmar Astrological Calendar Days

A typical Myanmar calendar usually describes astrological days such as auspicious days and ill fated days. In this post, we will present about the following days and javascript code for them.

CAN bus

CAN bus (controller area network) is a vehicle bus standard designed to allow microcontrollers and devices to communicate with each other. CAN bus is a message-based protocol, designed specifically for automotive applications but now also used in other areas such as aerospace, industrial automation and medical equipment. The advantages of CAN bus compared to RS232 communication are as follows.

Personal Release Form for Interviewer

Useful Web Resources for Academic Writing

Using Analog to Digital Converter of AT89C51CC01 Microcontroller

Using AT89STK-06 starter kit, I have written a few C code to read an analog to digital converter (ADC) input of AT89C51CC01 which is an 8051 microcontroller. It has 8 multiplexed ADC inputs with 10 bit resolution. As an example, ADC input pin 7 which is connected to a variable resister is read.

GNU Octave as an Alternative to MatLab

GNU Octave, FreeMat and Scilab are free open source software for numerical calculations alternative to MATLAB. After reading reviews in the Internet, I chose GNU Octave to try. Installing GNU Octave on a Windows system is easy. Download the files from its Download page and extract them to a folder. Then, you can put the shortcut at a desired location and change the path and the icon in its properties. When you open GNU Octave by clicking the shortcut, a window will appear as shown in the following figure. You can key in MATLAB commands there directly.

Chord-changer Javascript

I have written a Javascript program to transpose the guitar chords for a song into a different key. In the lyrics of the song, the guitar chords are supposed to be between ^and. You can try it at Myanmar Lyrics by clicking the Key Up and Key Down buttons. The codes are shown below.

Algorithm, Program and Calculation of Myanmar Calendar

မြန်မာဘာသာ ဖြင့်ဖတ်ရန်

An easier and faster modern method to calculate a past, present or future date in the Myanmar calendar (the Burmese calendar) is presented. The constants, formulas and steps are clearly defined to derive all the essential elements of the Myanmar calendar such as Myanmar year, Myanmar month, waxing or waning moon and day. Unlike the existing methods, this method is easy even for a person who is not familiar with the Myanmar calendar and its terms.

Non-inverting Amplifier Using op-amp

When PiezoDrive PDX 150 amplifier that I used to drive a piezoelectric actuator was damaged, I need to find a way to use an existing spare amplifier. The spare amplifier shown in the following figure has input voltage -2 V to 12 V and output voltage -20 V to 120 V. Its gain is only 10 while the damaged one is 20. Therefore, I decided to design and use an additional preamplifier with gain 2.


The input to that amplifier is -1 V to 6 V and its output shold be -2 V to 12 V. Since I need non-inverted output, I have used an op-amp with non-inverting configuration. I have also used a trimmer so that it is possible to fine tune the gain. Readily available LM358N op-amp IC is arbitrarily picked up. The schematic circuit diagram and resulting amplifier using available components in our lab is shown below. Since the non-inverting amplifier has very high input impedance, I have used additional 56k resister to prevent it from floating.


When tuning the amplifier to get exact gain of 2, it looked unreliable to use peak to peak readings from oscilloscope. That was why, I had to think of an alternative way to achieve it. I set channel 1 which was connected to output of the amplifier to 10 V per division and I set channel 2 which was connected to input of the implifier to 5 V per division. After that, I calibrated the trimmer to get identical waves at the screen :P

TRIAC Power Control Circuit

My vacuum cleaner was faulty. I thought the reason might be blown fuse and I tried to fix it. Unfortunately, I could not find any fuse either in machine or connector. I then thought that the reason might be the faulty power control circuit. Therefore, I tried to bypass the control circuit by directly connecting power and motor. But it did not work. Finally, I successfully solved the problem by throwing the old one and buying a new vacuum cleaner :P I kept the power control circuit and tried to trace it when I was free. It uses a simple RC circuit and a DIAC to make a phase lag in firing the TRIAC. The circuit is quite simple and it is shown below.

k-means clustering using custom distance measuring method

I have developed a MATLAB function to perform k-means clustering which enables custom distance measuring method. For example, to sort out histograms, chi-square distance may be more suitable. The following example function uses chi-square distance and you can replace it with any distance measurement method.

---

%k-means test program
X = [randn(100,2)+ones(100,2);...
     randn(100,2)-ones(100,2)];

[idx,ctrs] = KMeansCustom(X,2);
%[idx,ctrs] = kmeans(X,2);

plot(X(idx==1,1),X(idx==1,2),'r.','MarkerSize',12)
hold on
plot(X(idx==2,1),X(idx==2,2),'b.','MarkerSize',12)
plot(ctrs(:,1),ctrs(:,2),'kx',...
     'MarkerSize',12,'LineWidth',2)
legend('Cluster 1','Cluster 2','Centroids',...
       'Location','NW')

---

function [Idx,C]=KMeansCustom(X,k)
%KMeansCustom partitions the points in the n-by-d data matrix X into k clusters.
%[Idx,C]= KMeansCustom(X,k) returns 
%n-by-1 vector IDX containing the cluster indices of each point and 
%k-by-d matrix C containing the k cluster centroid locations.
%For n sample points with d dimensions in each point, X has n rows and d columns.
%File name: KMeanCustom.m
%Author: Yan Naing Aye
%Website: http://cool-emerald.blogspot.sg/


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Define maximum number of iterations
MaxIter=500;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[n,d]=size(X);
k=round(k);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%step1 :arbitrarily choose k samples as the initial cluster centers
p=randperm(n);
Mu=X(p(1:k),:);
D=zeros(k,d);
for t=1:MaxIter
 %step2:distribute the samples X  to the clusters 
 for j=1:k
        for i=1:n
            D(j,i)=ChiDist(X(i,:),Mu(j,:));%Use custom distance
        end
 end
 [ValMin,IndexMin]=min(D);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %step 3: update the cluster centers
    OldMu=Mu;
 for i=1:k
        Mu(i,:)=mean(X(IndexMin==i,:));
 end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %step4 :check convergence
 if sum(sum(abs(OldMu-Mu))) == 0 %< 1e-9
        break
 end
end
Idx=IndexMin';
C=Mu;

---

function d=ChiDist(v1,v2)
    dv=(v1-v2).^2;
    sv=abs(v1)+abs(v2);
    %------------------------------------------------------
    %eliminate zero denominator
    sv(sv==0)=1e-9;
    %------------------------------------------------------
    d=sum(dv./sv)./2;    
end

Amplitude spectrum of a signal using Fourier transform

I have developed a MATLAB function to calculate amplitude spectrum of a signal using Fourier transform. Example usage of the function and implementation is shown below. Another example MATLAB program using built in function fft is also demonstrated. I hope they will be helpful for those who want to find amplitude spectrum or power spectrum of a signal.

---

%Example MATLAB code to test function FX
Fs=100;
t=(0:1/Fs:1-1/Fs)';
x=2*cos(2*pi*5*t)+sin(2*pi*10*t);
[f a]=FX(x,Fs);
plot(f,a);

---

function [f a]=FX(x,fs)
%Calculates amplitude spectrum of a signal x using Fourier transform
%[f a]=FX(x,fs)
%Input: x=signal, fs = sampling rate
%Output: f = frequency axis, a = amplitude spectrum
%File name: FX.m
%Author: Yan Naing Aye
%Website: http://cool-emerald.blogspot.sg/

dT = 1/fs; 
N=length(x);
dF=1/(N*dT);
NF=fs/2;%Nyquist Freq
%You can always limit freq range for faster performance
%NF=20;
t=(0:dT:(N-1)*dT)';
f=(0:dF:NF)';
a=f;%initialize a
a(1)=mean(x);
for i=2 : length(a)
    b=(2*mean(x.*cos(2*pi*dF*(i-1)*t)));
    c=(2*mean(x.*sin(2*pi*dF*(i-1)*t)));
    a(i)=sqrt(b^2+c^2);
end

---

%MATLAB program to calculates amplitude spectrum of a signal x using fft
%Author: Yan Naing Aye
%Website: http://cool-emerald.blogspot.sg/

Fs=100;
t=(0:1/Fs:1-1/Fs)';
x=2*cos(2*pi*5*t)+sin(2*pi*10*t);
L=length(x);

A=2*abs(fft(x)/L);
A=A(1:Fs/2+1);
f = Fs/2*linspace(0,1,Fs/2+1);
plot(f,A);

---

Changing sampling rate using quadratic regression

I want to change sampling rate of a signal from 333Hz to 5000Hz. To get a smoother result, I have implemented second order hold system instead of popular zero order hold (ZOH) system. This approach can also be used for down sampling. Another possible application is to filter out noise without introducing phase delay. And it can be useful for real-time application.

Approach 1: Quadratic regression
Quadratic function is defined as

$$f=w_0+w_1 x + w_2 x^2$$

Then, cost function to be minimized is defined as

$$J(w)=\frac{1}{2}\sum_{i=1}^{n}(y_i-f_i)^2$$
$$J(w)=\frac{1}{2}\sum_{i=1}^{n}(y_i-w_0+w_1 x_i + w_2 x_i^2)^2$$

Optimal weights can be found by differentiating the cost function and setting them to zero.

$$\frac{\partial J(w)}{\partial w_0}=0$$
$$-\sum_{i=1}^{n}(y_i-w_0+w_1 x_i + w_2 x_i^2)=0$$
$$w_0\sum_{i=1}^{n}1+w_1\sum_{i=1}^{n}x_i+w_2\sum_{i=1}^{n}x_i^2=\sum_{i=1}^{n}y_i$$

Similarly, differentiating with w₁ and w₂ gives,

$$w_0\sum_{i=1}^{n}x_i+w_1\sum_{i=1}^{n}x_i^2+w_2\sum_{i=1}^{n}x_i^3=\sum_{i=1}^{n}y_i x_i$$
$$w_0\sum_{i=1}^{n}x_i^2+w_1\sum_{i=1}^{n}x_i^3+w_2\sum_{i=1}^{n}x_i^4=\sum_{i=1}^{n}y_i x_i^2$$

When these three equations are written in matrix form.

$$ \begin{bmatrix} \sum_{i=1}^{n}1 & \sum_{i=1}^{n}x_i & \sum_{i=1}^{n}x_i^2 \\ \sum_{i=1}^{n}x_i & \sum_{i=1}^{n}x_i^2 & \sum_{i=1}^{n}x_i^3 \\ \sum_{i=1}^{n}x_i^2 & \sum_{i=1}^{n}x_i^3 & \sum_{i=1}^{n}x_i^4 \end{bmatrix} \begin{bmatrix} w_0 \\ w_1 \\ w_2 \end{bmatrix} = \begin{bmatrix} \sum_{i=1}^{n}y_i \\ \sum_{i=1}^{n}y_i x_i \\ \sum_{i=1}^{n}y_i x_i^2 \end{bmatrix} $$
$$\mathbf{A}\mathbf{W}=\mathbf{B}$$
$$\mathbf{W}=\mathbf{A}^{-1}\mathbf{B}$$

Approach 2: Three equations
In our case, we can use only the last three points and we can get three equations from these three points.

$$y_1=w_0+w_1 x_1 + w_2 x_1^2$$
$$y_2=w_0+w_1 x_2 + w_2 x_2^2$$
$$y_3=w_0+w_1 x_3 + w_2 x_3^2$$
$$ \begin{bmatrix} 1 & x_1 & x_1^2 \\ 1 & x_2 & x_2^2 \\ 1 & x_3 & x_3^2 \end{bmatrix} \begin{bmatrix} w_0 \\ w_1 \\ w_2 \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} $$
$$\mathbf{A}\mathbf{W}=\mathbf{B}$$
$$\mathbf{W}=\mathbf{A}^{-1}\mathbf{B}$$

Approach 3: Two equations
If we define x₁=-1, x₂=0, and x₃=1, y₂ is equal to w₀. And,

$$y_1=y_2-w_1+ w_2$$
$$y_3=y_2+w_1+ w_2$$
$$ \begin{bmatrix} -1 & 1 \\ 1 & 1 \end{bmatrix} \begin{bmatrix} w_1 \\ w_2 \end{bmatrix} = \begin{bmatrix} y_1-y_2 \\ y_3-y2 \end{bmatrix} $$
$$ \begin{bmatrix} w_1 \\ w_2 \end{bmatrix} = \begin{bmatrix} -0.5 & 0.5 \\ 0.5 & 0.5 \end{bmatrix} \begin{bmatrix} y_1-y_2 \\ y_3-y2 \end{bmatrix} $$

Then, we have

$$w_0=y_2$$
$$w_1=-0.5(y_1-y_2)+0.5(y_3-y2)$$
$$w_2=0.5(y_1-y_2)+0.5(y_3-y2)$$

I have tested the above three approaches using MatLab and it is shown below.

%-------------------------------------------------------------------------
clc;
close all;
clear all;
%-------------------------------------------------------------------------
% y= w0 + w1*x + w2* x^2;
%-------------------------------------------------------------------------
%Got x and y
x=[-1 0 1]';
Wo=[4 3 2]';
y=Wo(1)+Wo(2)*x+Wo(3).*x.*x;

%-------------------------------------------------------------------------
%Approach 1
%Polynomial regression of order 2
%For n=3
S1=3;
Sx=sum(x);
Sx2=sum(x.*x);
Sx3=sum(x.*x.*x);
Sx4=sum(x.*x.*x.*x);
Sy=sum(y);
Syx=sum(y.*x);
Syx2=sum(y.*x.*x);

P=[S1 Sx Sx2; Sx Sx2 Sx3; Sx2 Sx3 Sx4];
B=[Sy Syx Syx2]';
%P1=P^(-1);
W1=P\B

%-------------------------------------------------------------------------
%Approach 2
%Linear equations
A=[1 x(1) x(1)*x(1);1 x(2) x(2)*x(2); 1 x(3) x(3)*x(3)];
W2=A\y
%-------------------------------------------------------------------------
%Approach 3
%Only 2 linear equations
w0=y(2);
w1=-0.5*( y(1)- y(2))+0.5*( y(3)- y(2));
w2=0.5*( y(1)- y(2))+0.5*(y(3)- y(2));
W3=[w0 w1 w2]'
%-------------------------------------------------------------------------

The following figure shows the result of using this method (blue color plot) compare to ordinary zero order hold (black color plot). This method gives much more smoother result but it should be noted that it introduces one sample delay.

The implementation of this method in LabVIEW using C code is shown in the following figure.

The first two approaches involve finding inverse of a 3x3 matrix and I have developed a C program as shown below.

#include
#include
main()
{
    float M[3][3]={{3,0,2},{0,2,0},{2,0,2}}; //initialize a 3x3 matrix
    float N[3][3]={{0,0,0},{0,0,0},{0,0,0}}; //allocate for inverse
    int i,j;
    float d;
    //-------------------------------------------------------------------------
    N[0][0]=(M[1][1]*M[2][2]-M[2][1]*M[1][2]);
    N[1][0]=-(M[1][0]*M[2][2]-M[2][0]*M[1][2]);
    N[2][0]=(M[1][0]*M[2][1]-M[1][1]*M[2][0]);
    d=M[0][0]*N[0][0]+M[0][1]*N[1][0]+M[0][2]*N[2][0];
    N[0][0]/=d;
    N[1][0]/=d;
    N[2][0]/=d;
    N[0][1]=-(M[0][1]*M[2][2]-M[0][2]*M[2][1])/d;
    N[1][1]=(M[0][0]*M[2][2]-M[0][2]*M[2][0])/d;
    N[2][1]=-(M[0][0]*M[2][1]-M[0][1]*M[2][0])/d;
    N[0][2]=(M[0][1]*M[1][2]-M[0][2]*M[1][1])/d;
    N[1][2]=-(M[0][0]*M[1][2]-M[0][2]*M[1][0])/d;
    N[2][2]=(M[0][0]*M[1][1]-M[0][1]*M[1][0])/d;
    //-------------------------------------------------------------------------
    //print 3x3 matrix
    for(i=0;i<3;i++)
    {
        for(j=0;j<3;j++) printf("%3.4f ",N[i][j]);
        printf("\n");
    }
    getch();
    return 0;
}

Rotations in 3D space using Euler angles

Rotations in 3D space to relate body reference frame of a device to world reference frame using an Euler angle sequence are discussed. Euler stated that

Any two independent orthonormal coordinate frames can be related by a sequence of rotations (not more than three) about coordinate axes, where no two successive rotations may be about the same axis.
The illustration of the world reference frame, {M}, and the body reference frame, {B}, is shown in the following figure.

Pan angle rotation, α, is acquired using computer vision. And, tilt angle rotation, β, and roll angle rotation, γ, are acquired by sensing the gravity, g, using onboard accelerometers. Let us consider the roll angle first, assuming there are no pan and tilt angles. When the roll angle, γ, is zero Y-axis is upward vertical and X-axis is horizontal pointing according to right handed rule. The gravity sensed by X-axis accelerometer and Y-axis accelerometer are denoted by gx and gy respectively. The definition of the roll angle is illustrated in the following figure.

The roll angle, γ, is calculated as

$$ \gamma = \tan^{-1}\frac{gx}{gy} $$.

Then, it is checked for second and third quadrants as

if(gy<0)
  if(gx>=0)
    γ=π+γ;
  else
    γ=-π+γ;
  end
end

The roll angle rotation matrix, R_γ, is obtained as follows.

$$ \mathbf{R}_{\gamma}= \begin{bmatrix} \cos(\gamma) & -\sin(\gamma) & 0 \\ \sin(\gamma) & \cos(\gamma) & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

After the roll angle is corrected, let us consider for tilt angle. The definition of the tilt angle is illustrated in the following figure.

The tilt angle, β, can be calculated from Y and Z components of the gravity. It is important to note that these gravity component values should be in the roll angle rotated frame. The new values for the gravity components are obtained as

$$ \begin{bmatrix} gx2 \\ gy2 \\ gz2 \end{bmatrix} = \mathbf{R}_{\gamma} \begin{bmatrix} gx \\ gy \\ gz \end{bmatrix} $$

The tilt angle, β, is calculated as
$$\beta = \tan^{-1}\frac{gy2}{gz2} $$.
Then, it is checked for second and third quadrants as

if(gz2>=0)
  if(gy2>=0)
    β=-π+β;
  else
    β=π+β;
  end
end

The tilt angle rotation matrix, R_β, is obtained as follows.

$$ \mathbf{R}_{\beta}= \begin{bmatrix} 1 & 0 & 0\\ 0 & \cos(\gamma) & -\sin(\gamma) \\ 0 & \sin(\gamma) & \cos(\gamma) \end{bmatrix}$$

Similarly, calculation of the pan angle rotation matrix, R_α, and description of pan angle definition are as follows.

$$ \mathbf{R}_{\alpha}= \begin{bmatrix} \cos(\alpha) & -\sin(\alpha) & 0 \\ \sin(\alpha) & \cos(\alpha) & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

The product of rotation matrices is itself a rotation matrix. $$ \mathbf{R}=\mathbf{R}_{\alpha} \mathbf{R}_{\beta} \mathbf{R}_{\gamma}$$
And, M=RB
Since R is a rotation matrix, inverse of R is obtained by transposing R. The implementation of this 3D rotation in LabVIEW using MabLab code is shown in the following figure.

Quaternion algebra is also easy and popular approach for such 3D transformation.

Reference:
Kuipers, Jack B., Quaternions and rotation sequences : a primer with applications to orbits, aerospace, and virtual reality, Princeton University Press, 1999, ISBN: 0691058725.