မြန်မာ ပြက္ခဒိန် အခါပေးရက်များ

မြန်မာ ပြက္ခဒိန်တွေမှာ ဗေဒင်နဲ့ ဆက်စပ်တဲ့ အခါပေးရက် ဥပမာ ရက်ကောင်း၊ရက်ဆိုး တွေကို ဖော်ပြလေ့ ရှိပါတယ်။ ဒီစာတမ်းမှာ အောက်ကရက်တွေနဲ့ သူတို့ကို javascript ကုဒ် သုံးပြီး ဘယ်လိုတွက်မလဲဆိုတာ ဖော်ပြပါမယ်

CAN bus

CAN Bus (controller area network) ဆိုသည်မှာ ဆက်သွယ်ရေး စနစ် တစ်ခုဖြစ်ပြီး မိုက်ခရိုကွန်ထရိုလာများ နှင့် အခြား အီလက်ထရွန်းနစ် ကိရိယာများကို ကွန်ပျူတာ မလိုအပ်ဘဲ အချင်းချင်း ဆယ်သွယ် နိုင်ရန် ပြုလုပ်ထားသော စနစ်ဖြစ်ပါတယ်။ CAN Bus ဟာ မက်ဆေ့ (message) အပေါ်တွင် အခြေခံသော ပရိုတိုကောလ်ဖြစ်ပြီး မော်တော်ယာဉ်များတွင် အဓိက အသုံးပြုလေ့ ရှိပေမယ့် စက်ရုံများနှင့် ဆေးဘက်ဆိုင်ရာ ကိရိယာများတွင်လည်း အသုံးပြုလေ့ ရှိပါတယ်။

Using Analog to Digital Converter of AT89C51CC01 Microcontroller

AT89STK-06 starter kit ကိုသုံးပြီး 8051 microcontroller တစ်ခုဖြစ်တဲ့ AT89C51CC01 အတွက် analog to digital converter (ADC) အဝင် ကိုဖတ်တဲ့ C ကုဒ် တစ်ချို့ကို ရေးထားပါတယ်။ အဲဒီ မိုက်ခရို ကွန်ထရိုလာ မှာ 10 bit resolution ရှိတဲ့ multiplexed လုပ်ထားတဲ့ အဝင် ၈ ခုရှိပါတယ်။ ဒီနမူနာမှာ variable resister နဲ့ ဆက်ထားတဲ့ အဝင် ပင်နံပါတ် ၇ ကို ဖတ်ကြည့်ထား ပါတယ်။

GNU Octave as an Alternative to MatLab

GNU Octave၊ FreeMat နဲ့ Scilab တို့က MATLAB လိုမျိုး တွက်ချက်မှုတွေ လုပ်ပေးနိုင်တဲ့ အလကားရတဲ့ free open source software တွေဖြစ်ပါတယ်။ လိုက်ရှာ ဖတ်ကြည့်တော့ GNU Octave က အကောင်းဆုံးလို့ ကြားတာနဲ့ သူ့ကို စက်ထဲမှာ ထည့်ပြီး စမ်းကြည့်တော့ တော်တော် သဘောကျ သွားပါတယ်။ GNU Octave ကို ဝင်းဒိုး ပေါ်မှာတင်ရတာ တော့ တော်တော်လွယ်ပါတယ်။ သူ့ရဲ့ Download page က ဖိုင်တွေကို ယူပြီး စက်ထဲမှာ folder တစ်ခု အနေနဲ့ ထုတ်ယူလိုက်ပါမယ်။ ပြီးရင် shortcut ဖိုင်ကို ကြိုက်တဲ့နေရာမှာ ထားပြီး သူ့ရဲ့ properties မှာ path နဲ့ icon ကို ပြင်ပေး လိုက်ရုံပါပဲ။ သူရဲ့ shortcut ကနေ Octave ကို ဖွင့်လိုက်ရင် အောက်က ပုံမှာ ပြထားတဲ့ အတိုင်း ကွန်မန်း ပေးလို့ရတဲ့ ဝင်းဒိုး ပွင့်လာပါမယ်။ အဲဒီမှာ MATLAB ကွန်မန်း တွေ တိုက်ရိုက် ရိုက်ထည့်ရင် လည်းရပါတယ်။

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

Read this article in English

မြန်မာ ပြက္ခဒိန် တွက်တဲ့ အခါ ပိုမို မြန်ဆန်လွယ်ကူ စေမယ့် ညီမျှခြင်းများကို တင်ပြမှာဖြစ်ပြီး၊ ပြီးခဲ့တဲ့ မြန်မာ ပြက္ခဒိန် ရက်ကိုပဲ ဖြစ်ဖြစ်၊ နောက်လာမည့် မြန်မာ ပြက္ခဒိန်ရက် ကို ကြိုပြီးပဲဖြစ်ဖြစ် ဘယ်လို တွက်မလဲ ဆိုတာ ဆွေးနွေးပါမယ်။ မြန်မာ ပြက္ခဒိန် ရက်စွဲ တစ်ခုရဲ့ မြန်မာ ခုနှစ်၊ မြန်မာလ၊ လဆန်း လဆုတ်၊ မြန်မာရက်၊ အဲဒီနှစ်က ဝါထပ် မထပ်၊ ထပ်ရင် လည်း ဝါကြီးလား၊ ဝါငယ်လား ဆိုတာ ကိုအလွယ်တကူ တွက်ထုတ်နိုင်ဖို့ ကိန်းသေတွေ၊ ဖော်မြူလာတွေ၊ တွက်ချက်ပုံ အဆင့် တွေ ကို ရှင်းရှင်းလင်းလင်း တင်ပြပါမယ်။ မြန်မာ ပြက္ခဒိန်၊ နတ္ခတ် အခေါ်အဝေါ် တွေနဲ့ မရင်းနှီးရင် တောင်မှ အလွယ် တစ်ကူ နားလည် နိုင်မှာပါ။

Non-inverting Amplifier Using op-amp

လတ်တလော သုံးနေတဲ့ amplifier (PiezoDrive PDX 150) လေးက ပျက်သွားတော့ အပိုရှိနေတဲ့ စပယ်ယာ amplifier လေးနဲ့ ပြောင်းသုံးဖို့ လိုလာပါတယ်။ အောက်ကပုံမှာ ပြထားတဲ့ စပယ်ယာ amplifier လေးက အဝင်ဗို့ -2 V ကနေ 12 V ရှိပြီး၊ အထွက်ဗို့က -20 V ကနေ 120 V ထုတ်ပေးနိုင်ပါတယ်။ ပျက်သွားတဲ့ ကောင်က ချဲ့ဆ ၂၀ ရှိပြီး၊ အပို ရှိတဲ့ amplifier လေးရဲ့ ချဲ့ဆ က ၁၀ ဆပဲရှိပါတယ်။ အဲဒီတော့ ချဲ့ဆ ၂ ဆရှိတဲ့ preamplifier လေး တစ်ခု ကို ကိုယ့်ဟာ ကိုယ် ဒီဇိုင်းလုပ်ပြီး ရှေ့ကနေ ခံသုံးဖို့ ဆုံးဖြတ်လိုက်ပါတယ်။

လုပ်ရမယ့် amplifier ရဲ့ အဝင်ဗို့က -1 V ကနေ 6 V ရှိမှာ ဖြစ်ပြီး၊ အထွက်ဗို့က -2 V ကနေ 12 V ထုတ်ပေးဖို့ လိုပါမယ်။ အဝင်ဗို့ နဲ့ အထွက်ဗို့က အတည့် (non-inverted) ဖြစ်ဖို့ လိုတာမို့ op-amp ကိုnon-inverting configuration သုံးပြီး ဒီဇိုင်းလုပ်လိုက်ပါတယ်။ ချဲ့ဆကို အသေးစိတ် ချိန်ညှိ ချင်တာဖို့ trimmer လေးကို ပါထည့်သုံးလိုက်ပါတယ်။ အလွယ်တကူ တွေ့တဲ့ LM358N op-amp IC ကိုပဲ ကောက်သုံးလိုက်ပါတယ်။ ဆားကစ် ပုံ နဲ့ lab ထဲမှာ အလွယ်တကူ ရတဲ့ ပစ္စည်းလေးတွေသုံးပြီး ရလာတဲ့ ဟာလေးကို အောက်မှာ ပြထားပါတယ်။ ဒီဆားကစ် ရဲ့ အဝင် ခုခံမှု (input impedance) က အရမ်းကြီးတာမို့ အဝင်ကို ဘာမှ ဆက်မထားတဲ့ အချိန် မှာ floating ဖြစ်မနေအောင် 56k ရီစစ္စတာလေး ထည့်လိုက်ပါတယ်။


ဒီ amplifier ရဲ့ ချဲ့ဆ ကို ၂ ဆတိတိရအောင် ချိန်တဲ့ အခါ oscilloscope ရဲ့peak to peak တန်ဖိုး ပြနေတာ ကို ဖတ်ရတာ သိပ်မသေချာဘူး ထင်တာနဲ့ တစ်ခြား နည်းလမ်း ကို ထပ်စဉ်းစား ကြည့်ပါတယ်။ နောက်တော့ ချန်နယ် ၁ နဲ့ ဆက်ထားတဲ့ အထွက် ကို 10 V per division ထားပြီး၊ ချန်နယ် ၂ နဲ့ ဆက်ထားတဲ့ အဝင်ကို 5 V per division ထားလိုက်ပါတယ်။ ပြီးတာနဲ့ trimmer လေးကို ဖြည်းဖြည်းခြင်း လှည့်ပြီး ထပ်တူကျတဲ့ လှိုင်းနှစ်ခု ရအောင် ချိန်လိုက်ပါတယ်။

TRIAC Power Control Circuit

အိမ်က ဖုန်စုပ်စက်လေး ပျက်သွားတော့ ဖျု့စ် ပျက်တာလောက် ဆိုရင် ပြင်လို့ ရရင် ရနိုင်ကောင်းရဲ့ ဆိုပြီး စက်ကို ဖွင့် ဖျု့စ်လိုက်ရှာတော့ မတွေ့ပါဘူး :P။ ပလပ်မှာ လည်း မပါတော့ ဒါ ဆိုရင် ပါဝါ ကွန်ထရိုး လုပ်ပေးတဲ့ ဆားကစ် ကပျက်တာ လည်း ဖြစ်နိုင်တာပဲ ဆို ပြီး အဝင်ကြိုးကို ဆားကစ် ကို မဆက်ပဲ ကျော်ပြီး မော်တာ ကို တိုက်ရိုက် ဆက်ကြည့်ပါတယ်။ ဒါလည်း မလုပ်ပါဘူး။ နောက်ဆုံး တော့ စက် အဟောင်းကို ပစ်ပြီး၊ ဖုန်စက် စက် အသစ် တစ်လုံး ဝယ် ပြီး ဖြေရှင်းလိုက်တယ်။ အေးရော :P။ ပါဝါ ကွန်ထရိုး လုပ်ပေးတဲ့ ဆားကစ် လေးကိုတော့ မပစ်ပဲ အားတဲ့အချိန် trace လိုက်ကြည့်ပါတယ်။ သူက ရိုးရှင်းတဲ့ RC ဆားကစ်လေးနဲ့ DIAC ကိုသုံးပြီး TRIAC ဘယ်လောက် နောက်ကျပြီးမှ ပွင့် မလဲဆိုတဲ့ အချိန်ကို ထိန်းသွားတာပါ။ ဆားကစ်လေးက ရိုးရှင်းပြီး အခြေခံ ကျပါတယ်။ သူ့ကို အောက်ကပုံတွေမှာ ပြထားပါတယ်။

k-means clustering using custom distance measuring method

ကြိုက်တဲ့ အကွာအဝေး တိုင်းတဲ့ နည်းသုံးပြီး k-means clustering လုပ်ပေးတဲ့ MATLAB ဖန်ရှင် လေး တစ်ခု ရေးကြည့် ထားပါတယ်။ ဥပမာ histograms တွေကို ခွဲခြားဖို့ဆိုရင် chi-square distance ကို သုံးတာက ပိုကောင်း ချင် ကောင်းနိုင်ပါတယ်။ အောက်က နမူနာမှာ chi-square distance သုံးထားတာကို ပြထားပေမယ့် အဲဒီ distance ဖန်ရှင်နေရာမှာ ကိုယ်ကြိုက်တဲ့ နည်း ကို အစားထိုးပြီး သုံးလို့ရပါတယ်။

---

%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

signal တစ်ခုရဲ့ amplitude spectrum ကို Fourier transform သုံးပြီး ရှာတဲ့ MATLAB ဖန်ရှင်လေး တစ်ခု ရေးကြည့်ထားပါတယ်။ အဲဒီ ဖန်ရှင်ကို နမူနာ သုံးတဲ့ ပုံစံ ရယ်၊ ဖန်ရှင် ရေးထားတဲ့ ကုဒ်တွေ ရယ်ကို အောက်မှာ ဖော်ပြထားပါတယ်။ MATLAB မှာ ပါတဲ့ ဖန်ရှင် fft ကို သုံးတဲ့ ဟာ ကို လည်း နောက်ထပ် နမူနာ တစ်ခု ပြထားပါတယ်။ ဒီဟာလေးတွေက amplitude spectrum တို့၊ power spectrum တို့ ရှာချင်တဲ့ သူတွေ အတွက် အသုံးဝင်မယ်လို့ မျှော်လင့်ပါတယ်။

---

%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

signal တစ်ခုရဲ့ sampling rate ကို 333Hz ကနေ 5000Hz ကို ပြောင်းဖို့ အတွက် လုပ်ကြည့်ထားပါတယ်။ ပိုမို ချောမွေ့တဲ့ ရလာဒ် ရချင်တာမို့ ပုံမှန် zero order hold (ZOH) ကို သုံးမယ့်အစား second order hold ကို သုံးကြည့်ထားတာပါ။ ဒီမှာ sampling rate ကို တင်ဖို့လုပ်ထားပေမယ့်၊ ဒီနည်းကို down sampling လုပ်ဖို့လည်း သုံးလို့ရပါတယ်။ အခြား ဖြစ်နိုင်ချေတွေကတော့ noise ကို phase delay မရှိပဲ filter လုပ်တဲ့ဟာမျိုးတွေ အတွက်လည်း သုံးလို့ ရနိုင်ပါတယ်။ real-time application တွေအတွက်ပေါ့။
Approach 1: Quadratic regression
Quadratic ဖန်ရှင်တစ်ခုကို အောက်ပါအတိုင်း သတ်မှတ်နိုင်ပါတယ်။

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

အနည်းဆုံး ဖြစ်အောင် လုပ်မယ့် cost function ကို အောက်ကလို သတ်မှတ်မှာပါ။

$$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 တွေကို differentiate လုပ်ပြီးသုည ညီလိုက်ပြီး ရှာနိုင်ပါတယ်။

$$\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$$

w₁ နဲ့ w₂ ကိုလည်း အဲဒီလိုပဲ ရှာလိုက်မယ် ဆိုရင် အောက်အတိုင်း ထပ်ရလာပါမယ်။

$$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$$

အဲဒီ equation သုံးခုကို မေ့ထရစ် ပုံစံ နဲ့ အောက်ပါအတိုင်း ရေးလို့ရပါတယ်။

$$ \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
တကယ်တော့ ဒီကိစ္စမှာ နောက်ဆုံး အမှတ် ၃ ခုကိုပဲ သုံးလိုက်ရင်ရတဲ့ အတွက်၊ ညီမျှခြင်း သုံးကြောင်း ကို အောက်က အတိုင်း တစ်ခါတည်း ချရေးလို့လည်း ရပါတယ်။

$$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 အကယ်၍ x₁=-1, x₂=0, နဲ့ x₃=1 လို့ သတ်မှတ်လို့ ရတယ်ဆိုရင် y₂ က w₀ နဲ့ ညီသွားပါမယ်။ နောက်ကျန်တဲ့ အမှတ် ၂ ခုကနေ ညီမျှခြင်း နှစ်ကြောင်းပဲ ရှင်းဖို့ လိုပါတော့တယ်။

$$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} $$

နောက် အောက်ကအတိုင်းရလာပါမယ်။

$$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)$$

အထက်ကနည်းတွေကို MatLab နဲ့စမ်းကြည့်ထားတာကို အောက်က ကုဒ်တွေမှာ တွေ့နိုင်ပါတယ်။

%-------------------------------------------------------------------------
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]'
%-------------------------------------------------------------------------

အောက်ကပုံမှာ ဒီနည်းကို သုံးလို့ရတဲ့ ရလဒ် (အပြာ) နဲ့ ပုံမှန် zero order hold သုံးရင် ရမယ့် ရလဒ် (အနက်) တို့ကို ယှဉ်ပြထားပါတယ်။ ဒီနည်းက ပိုပြီး ညက်ညော ချောမွေ့တဲ့ ရလဒ်ကို ရပေမယ့် one sample delay ရှိသွားတာကိုတော့ သတိထားဖို့ လိုပါတယ်။

ဒီနည်းလမ်းကို LabVIEW မှာ C code သုံးပြီး လုပ်ထားတာကို အောက်ကပုံမှာ တွေ့နိုင်ပါတယ်။

ပထမ နှစ်ခုမှာ 3x3 matrix ကို inverse ရှာဖို့ လိုတဲ့အတွက် inverse ရှာတဲ့ C program လေးတစ်ခုကို ရေးပြီး စမ်းကြည့်ထားပါတယ်။

#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

ကိရိယာ တစ်ခုရဲ့ body reference frame နဲ့ world reference frame ရဲ့ ဆက်သွယ်မှု ကို ရှာဖို့အတွက် Euler angle sequence တစ်ခုကို သုံးပြီး သရီးဒီ စပေ့ ထဲမှာ rotation လုပ်တဲ့ အကြောင်း ပြောချင်ပါတယ်။ အွိုင်လာ တင်ပြခဲ့တဲ့ သီအိုရမ်တစ်ခု မှာ ပြောထားတာက အောက်ပါအတိုင်း ဖြစ်ပါတယ်။
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.
{M} နဲ့ ကိုယ်စားပြု ခေါ်မယ့် world reference frame နဲ့ {B} နဲ့ ကိုယ်စားပြုမယ့် body reference frame တို့ကို အောက်မှာ ပုံနဲ့ နမူနာ ပြထားပါတယ်။

Pan angle rotation ကို α လို့ ခေါ်မှာ ဖြစ်ပြီး သူ့ကို computer vision ကနေ ရနိုင်ပါတယ်။ နောက် tilt angle rotation - β နဲ့ roll angle rotation - γ တို့ကိုတော့ ကိရိယာ ထဲမှာ တပ်ထားမယ့် accelerometers တွေကနေ မြေဆွဲအား (gravity- g ) ကို အာရုံခံ ပြီး သိရှိနိုင်ပါတယ်။ အရင်ဆုံး pan and tilt မရှိဘူး လို့ ယာယီ ယူဆပြီး roll angle ကို စဉ်းစားကြည့်ပါမယ်။ γ က သုည ဖြစ်နေချိန်မှာ Y ဝင်ရိုးက အထက်ကို ညွှန်ပြီး၊ X ဝင်ရိုးက right handed rule အတိုင်း ရေပြင်ညီ ရှိမယ်လို့ သတ်မှတ် ပါမယ်။ X ဝင်ရိုး accelerometer နဲ့ Y ဝင်ရိုး accelerometer တို့က မြေဆွဲအားကို အာရုံခံ တိုင်းလို့ရတဲ့ တန်ဖိုးတွေကို gx နဲ့ gy လို့ခေါ်ပါမယ်။ roll angle ကို သတ်မှတ် ထားပုံကို အောက် က ပုံမှာ ဖော်ပြထားပါတယ်။

Roll angle - γ ကို အောက်က ညီမျှခြင်း ကိုသုံးပြီး တွက်ယူလို့ ရပါတယ်။

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

နောက် second နဲ့ third quadrants တွေ မှာ ရှိမရှိကို ပြန်စစ်ပါမယ်။

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

ဒါဆို roll angle rotation matrix - R_γ ကို အောက်က အတိုင်းရနိုင်ပါတယ်။

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

Roll angle ကို ပြင်ပြီးတဲ့နောက် tilt angle ကိုဆက်စဉ်းစားပါမယ်။ Tilt angle ရဲ့ သတ်မှတ်ပုံ ကို အောက်မှာ ပြထားပါတယ်။

Tilt angle - β ကို Y နဲ့ Z မြေဆွဲအား တိုင်းတာမှုတွေ ကနေ တွက်လို့ရပါတယ်။ သတိထား ဖို့ လိုတာကတော့ အဲဒီတန်ဖိုးတွေက roll angle ပြင်ထားတဲ့ တန်ဖိုးတွေ ဖြစ်ရပါမယ်။ ပြောင်းသွားတဲ့ Y နဲ့ Z မြေဆွဲအား တိုင်းတာမှုတွေ ကို အောက်ပါအတိုင်း တွက်ယူနိုင်ပါတယ်။

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

ဒါဆို β ကို အောက်ပါ ညီမျှခြင်း ကနေရှာလို့ ရပါပြီ။

$$ \beta = \tan^{-1}\frac{gy2}{gz2} $$.

နောက် second and third quadrants တွေ အတွက် စစ်ပါမယ်။

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

အဲဒီအခါ Tilt angle rotation matrix -R_β က အောက်ပါအတိုင်းရလာပါမယ်။

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

အဲဒီလိုပါပဲ၊ pan angle rotation matrix - R_α နဲ့ pan angle ရဲ့ သတ်မှတ်ပုံတွေ ကိုလည်း အောက်ပါအတိုင်း တွေနိုင်ပါတယ်။

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

အဲဒီ rotation matrices တွေရဲ့ မြှောက်လဒ် က rotation matrix တစ်ခု ပါပဲ။ $$ \mathbf{R}=\mathbf{R}_{\alpha} \mathbf{R}_{\beta} \mathbf{R}_{\gamma}$$
အဲဒီအခါ အောက်ပါအတိုင်း ရရှိလာမှာပါ။
M=RB

R က rotation matrix တစ်ခု ဖြစ်တာမို့ သူ့ရဲ့ inverse ကို ရှာဖို့ လိုခဲ့ရင် အဲဒါက R transpose နဲ့ အတူတူပါပဲ။ ဒီ 3D rotation ကို LabVIEW မှာ MabLab code သုံးပြီး လုပ်ထားတာကို အောက်ကပုံမှာ မြင်နိုင်ပါတယ်။

အခုလို 3D transformation တွေအတွက် Quaternion algebra ကလည်း လွယ်ကူပြီး ခေတ်စားနဲ့ နည်းတစ်နည်း ပါပဲ။

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