เริ่มการวิเคราะห์หุ้นด้วย MATLAB

เห็นช่วงนี้หุ้นผันผวนมาก เลยคิดว่าควรจะศึกษาแนวทาง technical ไว้ ตอนนี้ยังไม่มีความรู้อะไร แต่ก่อนจะวิเคราะห์ได้ เราต้องมีข้อมูลก่อน ทีนี้เราจะดึงข้อมูลหุ้นมายังไงล่ะ วันนี้จะลองมาทำให้ดู อาจจะไม่ใช่วิธีที่ดีที่สุด แต่มันเวิร์คก็ok

ก่อนอื่นเราต้องมี python 3.6+  มาก่อน เพราะว่ายังหาวิธีดึงข้อมูลด้วย matlab ไม่เป็น 555 หลังจากนั้นก็ลง pandas_datareader module ใน command prompt ด้วยคำสั่ง

pip install pandas-datareader

เสร็จแล้วก็เปิด python IDE ขึ้นมาแล้วพิมพ์คำสั่งเพื่อที่จะดึงข้อมูลหุ้นได้เลย สมมติว่าเราจะดึงข้อมูลของ PTTGC ณ วันนี้ 1 มกรา ถึง 24 กรกฎา 2020 แล้วก็exportข้อมูลไปเป็นไฟล์ csv ก็พิมพ์ตามนี้

import pandas_datareader as pdr
import pandas as pd
import datetime

pttgc = pdr.get_data_yahoo('PTTGC.BK',
                           start=datetime.datetime(2020, 1, 1),
                           end=datetime.datetime(2020, 7, 24))
df = pd.DataFrame(pttgc)
df.to_csv (r'E:\backupusb\myMATLAB\export_pttgc.csv', header=True)
print(df)

เมื่อทำการเซฟข้อมูลลงเครื่องเรียบร้อย เราก็มาทำการplotใน MATLAB ได้

ถามว่าทำไมถึงไม่ใช้pythonทำ...คำตอบง่ายมาก "ก็สะดวก MATLAB กว่า" นั่นเอง 5555

เอ้า เปิด MATLAB ขึ้นมาแล้วพิมพ์คำสั่งไปตามนี้

clear;
filename = 'export_pttgc.csv';
[num,txt,raw] = xlsread(filename);
%Date|High|Low|Open|Close|Volume|Adjust close
hi = num(:,1);
lo = num(:,2);
op = num(:,3);
cl = num(:,4);
vol = num(:,5);
adjcl = num(:,6);
time = datetime(txt(2:end,1),'InputFormat','MM/dd/yyyy');
plot(time,hi,'-b','LineWidth',2); hold on;
plot(time,lo,'-r','LineWidth',2); legend('high','low');
axis tight;hold off;
title('pttgc')

เราก็จะได้กราฟสวยมาเริ่มเล่นใน episode ต่อไป ...

ตอนนี้ขอตัวไปเรียนการวิเคราะห์หุ้นแบบ technical ก่อน มีอะไรใหม่เดี๋ยวจะมาอัพเดทให้

สร้างจำนวนเต็มที่ไม่ซ้ำกัน

อันนี้ได้ไอเดียมาตอนเล่นเกมส์แล้วมีเพื่อนในกิลด์ถามว่า จะgenerateเลขจำนวนเต็มให้ไม่ซ้ำกันได้ยังไง พอดีเค้าเพิ่งหัดเรียนเขียนโปรแกรม ไอ้เราก็ไม่ได้พูดอะไรมาก แต่ก็กลับมานั่งคิด เอ ถ้าเป็นเราจะทำยังไงน้า เลยคิดชุดคำสั่งง่ายๆขึ้นมา มันคงไม่ตอบโจทย์ที่เค้าถาม แต่เห็นว่าน่าสนใจดีเลยเอามาแปะไว้

สมมติเราอยากได้จำนวนเต็มไม่ซ้ำกันภายใน 1 ถึง 1 แสน จำนวนหมื่นตัว ก็สร้างเลขจำนวนเต็มขึ้นมาก่อนด้วยคำสั่ง randi

nid = 1e4;
uid = randi([1,1e5],1,nid);

ทีนี้มันต้องมีเลขที่ซ้ำกันบ้างล่ะ เราก็จะมาคัดเลขออกกัน ถ้าเป็นคนอื่นจะทำยังไง ส่วนใหญ่คงใช้ for loop แต่ไม่เอา เราอยากได้วิธีอื่น ที่คิดขึ้นมาเอง 555

[val,ind] = sort(uid);
dupcheck = find(diff(val)==0);
val(dupcheck) = [];
ind(dupcheck) = [];

อันดับแรก เรียง vector จากจำนวนมากไปน้อยก่อน แล้วก็หาว่าตัวเลขที่อยู่ติดกันเหมือนกันหรือไม่ ด้วยการหา different ใน vector นี้ แล้วหาตำแหน่งที่มีค่าเท่ากับ 0 เพราะตำแหน่งที่มีค่าเท่ากับศูนย์คือตำแหน่งที่ตัวเลขเหมือนกัน หลังจากได้ตำแหน่ง dupcheck มาแล้วก็ไปตัดตัวเลขที่ซ้ำกันออก ใน val และ ind

newid = nan(1,nid);
newid(ind) = uid(ind); 
newid(isnan(newid))=[];
plot(newid,'ob');

หลังจากนั้นสร้างตัวแปรใหม่ newid ขึ้นมา โดยให้มีขนาดเท่ากับ nid แล้วก็ทำการยัดตัวเลขจาก nid เข้า newid โดยตัดตำแหน่งที่ซ้ำกันออก
เท่านี้เราก็จะได้จำนวนเต็มที่ไม่ซ้ำกันแล้ว เพียงแต่วิธีนี้จำกัดจำนวนตัวเลขได้เท่าที่ memory ของ computer จะเก็บ nid ไหว และที่สำคัญมันบอกจำนวนที่จะgenerateออกมาไม่ได้ เช่นถ้าบอกว่าอยากได้ unique number 9000 ตัว ไรงี้ ก็จะทำให้ไม่ได้ แต่ก็ ไม่เลวร้ายเกินไปสำหรับ quiz ลับสมองนี้ล่ะมั้ง

การหาจุดตัดระหว่างเส้นตรงและวงกลม MATLAB version

จริงๆแล้ว MATLAB มันมีชุดคำสั่งในการหาจุดตัดนี้อยู่นะ มันคือ linecirc ตัวอย่างเช่น

slope = 0;
intercpt = 0.3;
centerx = 0;
centery = 0;
radius = 0.5;
[xout,yout] = linecirc(slope,intercpt,centerx,centery,radius)

คำตอบที่ได้คือ
xout =

    0.4000   -0.4000


yout =

    0.3000    0.3000

ส่วนค่าในแต่ละช่องคืออะไร เราจะมารีวิวสมการเส้นตรงและสมการวงกลมกันวันนี้

ก่อนอื่น สมการวงกลม

\[ (x-u)^2 + (y-v)^2 = r^2\]

โดยที่ (u,v) คือตำแหน่งจุดศูนย์กลางวงกลม แทนด้วย centerx,centery ใน MATLAB code ข้างบน
r คือ รัศมีวงกลม แทนด้วย radius ในโค้ดด้านบน

ส่วนสมการเส้นตรงแทนด้วย

\[ y = mx+c\]

โดย m คือความชันของกราฟเส้นตรง และ c คือจุดตัดแกน y

ตำแหน่งที่เส้นตรงตัดกับวงกลม คือตำแหน่งที่ x,y ของสมการวงกลม มีค่าเท่ากับ x,y ของสมการเส้นตรง ดังนั้นเราจึงจัดเรียงสมการใหม่ โดยการแทนค่า y ของสมการเส้นตรงในสมการวงกลม จะได้สมการดังนี้
\[ (x-u)^2 + (mx+c-v)^2 = r^2\]
ทีนี้ก็จัดเรียงสมการใหม่เพื่อที่จะหาค่า x ก่อนอื่นกระจายตัวแปร polynomial ก่อน

\[ x^2 - 2xu + u^2 + (mx)^2 +2mx(c-v) + (c-v)^2 = r^2\] \[ x^2 - 2xu + u^2 + (mx)^2 +2mxc-2mxv + c^2 -2cv + v^2 = r^2\] เสร็จแล้วก็จัดกลุ่มใหม่เพื่อให้ดูง่าย \[ [x^2+ (mx)^2] +[-2xu+2mxc-2mxv] + [u^2 + c^2 -2cv + v^2 - r^2] = 0\] \[ [1+m]x^2 +[-2u+2mc-2mv]x + [u^2 + c^2 -2cv + v^2 - r^2] = 0\] หลังจากนั้นก็หาค่า x โดยใช้คำสั่ง roots ใน MATLAB จะได้ดังนี้

m = 0;
c = 0.3;
u = 0;
v = 0;
r = 0.5;

a1 = (1+m)^2;
a2 = -2*u + 2*m*c - 2*m*v;
a3 = u^2 + c^2 - 2*c*v + v^2 - r^2;
x = roots([a1,a2,a3])';

คำตอบ x ที่ออกมาคือ
x =

    0.4000   -0.4000

เหมือนกันกับข้างบนเลย ทีนี้พอเราได้ค่า x แล้วก็สบายล่ะ เอาค่าไปแทนในสมการเส้นตรงหาค่า y ได้เลย

y = m.*x+c

คำตอบที่ได้คือ
y =

    0.3000    0.3000
ที่วันนี้มาทำวิธีmanualให้ดูแทนที่จะใช้คำสั่ง linecirc เพราะว่าถ้าเรารู้วิธีการคิดเราจะสามารถนำไปต่อยอดได้ โดยเฉพาะคำสั่งของ MATLAB นี้ ไม่ครอบคลุมถึงรูป 3 มิติ เพราะงั้นรู้ที่มาไว้ไม่เสียหาย

การหา local maxima และ minima อย่างรวดเร็ว

สมมติว่าเรามีสัญญาณหนึ่ง แล้วเราต้องการหาว่าจุดยอดอยู่ที่ไหน คำสั่งใน MATLAB ที่ใช้เป็นปกติคือ findpeaks แต่ถ้าใครเคยใช้คำสั่งนี้ใน loop จะเห็นว่ามันช้ามาก ดังนั้นสำหรับฟังชั่นที่ไม่มีสัญญาณรบกวน มันมีวิธีอื่นที่ง่ายและเร็วกว่า วันนี้จะมาแสดงวิธีให้ดู

ก่อนอื่น สมมติเรามี sinc function signal ละกัน

x = linspace(-10,10,501);
y = sinc(x);

เราสามารถหาความชันของสัญญาณด้วยคำสั่ง diff

A = diff(y);

เมื่อเทียบกันสัญญาณรูปบน จะเห็นว่าจุดสูงสุดของสัญญาณคือจุดที่ค่าควาามชันเปลี่ยนจาก + เป็น - และจุดต่ำสุดของสัญญาณคือจุดที่ความชันเปลี่ยนจาก - เป็น + เพราะว่าเราสนใจแค่เครื่องหมาย เราจึงเปลี่ยนค่าความชันเป็นเครื่อง + หรือ - ด้วยคำสั่ง sign

A = sign(A)

หลังจากนั้นเราจะหาตำแหน่งที่สัญญาณเปลี่ยนเครื่องหมายจาก + เป็น -  และ - เป็น + ด้วยคำสั่ง conv

A = conv(A,[-1,1]);

ทีนี้เราก็จะหาจุดสูงสุดได้เมื่อ A = 2 และ จุดต่ำสุดที่ A = -2

maximas = find(A == 2);
minimas = find(A == -2);

เพื่อความกระชับ เราสามารถเขียนโค้ดทั้งหมดได้ในบรรทัดเดียวดังนี้

maximas = find(conv(sign(diff(y)),[-1,1]) == 2);
minimas = find(conv(sign(diff(y)),[-1,1]) == -2);

เพื่อยืนยันว่ามันใช้ได้เดี๋ยวจะขอพล๊อตกราฟให้ดู ดังนี้

x = linspace(-10,10,501);
y = sinc(x);
ymax = NaN(length(y));
ymin = ymax;
maximas = find(conv(sign(diff(y)),[-1,1]) == 2);
minimas = find(conv(sign(diff(y)),[-1,1]) == -2);
ymax(maximas) = y(maximas);
ymin(minimas) = y(minimas);
plot(x,y,'-b','LineWidth',2);hold on;
plot(x,ymax,'*r','LineWidth',2);
plot(x,ymin,'or','LineWidth',2);
set(gca,'FontSize',16);
hold off;

การคูณเมตริกซ์แบบต่างๆ

ในโปรแกรม MATLAB เป็นที่รู้กันว่ามีการคูณmatrix 2 แบบ

1. Element wise multiplication แทนด้วย A.*B คือการคูณแบบตำแหน่งต่อตำแหน่ง
2. Matrix multiplication แทนด้วย A*B คือการคูณเมตริกซ์

โพสที่แล้วเราโชว์การคูณแบบ element wise ให้ดู โพสนี้เราจะมาแสดงการคูณแบบเมตริกซ์ให้ดู 

แต่ก่อนอื่นเรามาทบทวนนิยามการคูณแบบเมตริกซ์ก่อน

สมมติเรามี matrix A และ B ดังนี้

ผลคูณ A*B = C จะได้

โดยที่ 

การคูณใน MATLAB อาจทำได้ 3 แบบ แต่วันนี้จะลองแค่ 2 แบบให้ดูเป็นตัวอย่าง

A = [2,3;3,5];
B = [3,5;2,7];
tic
C1 = A*B;
toc
tic
C2 = sum(bsxfun(@times, permute(A,[1,3,2]),permute(B,[3,2,1])),3);
toc

C1 ใช้เวลา 0.000134 seconds.
C2 ใช้เวลา 0.000244 seconds.

แบบแรกใช้คำสั่งคูณธรรมดา ใช้เวลาน้อยกว่าแบบที่2 ที่ให้ bsxfun .... ก็สรุปว่าการคูณแบบธรรมดา เร็วกว่าการคูณแบบพิสดารนะคะ

แต่ค่ะแต่...เราคงไม่แสดงวิธีที่สองให้ดูถ้ามันไม่มีอะไรน่าสนใจ
สำหรับ 2D matrix การคูณธรรมดานั้นเร็วที่สุด แต่ถ้าสมมติว่าเรามี 2D matrix หลายๆอันมาต่อกันในมิติที่สาม, k,แล้วเราต้องการคูณเมตริกซ์ k ครั้ง เราควรจะใช้วิธีไหน

โค้ดที่ใช้ในการทดสอบเป็นดังนี้

n = 3;
A = rand(n,n);
B = rand(n,n);
[m,n] = size(A);
vecK = [10,50,100,200,500,1000,3000,1e4];

t1 = zeros(1,length(vecK));
t2 = t1;
t3 = t1;
for ii = 1:length(vecK)
    K = vecK(ii);
    k = 1:K;
    A1 = ones(1,1,K).*A;
    B1 = ones(1,1,K).*B;
    C1 = zeros(m,n,K);
    tic;
    for k = 1:K
      C1(:,:,k) = A1(:,:,k) * B1(:,:,k);
    end
    t1(ii)=toc;

    tic;
    C2 = zeros(m,n,k);
    for i = 1:m
      for j = 1:m
        for k = 1:n
          C2(i,k,:) = C2(i,k,:) + A1(i,j,:).*B1(j,k,:);
        end
      end
    end
    t2(ii)=toc;

    tic;
    C3 = zeros(m,n,K);
    for j = 1:size(A,2)
      C3 = C3 + bsxfun(@times,A1(:,j,:),B1(j,:,:));
    end
    t3(ii)=toc;
end
h=plot(vecK,t1,'-*b',vecK,t2,'-*r',vecK,t3,'-*g');
h(1).LineWidth = 2;
h(2).LineWidth = 2;
h(3).LineWidth = 2;
set(gca,'FontSize',14);
axis tight;
xlabel('k');ylabel('time(s)');grid on;
legend('A*B', 'vectorized', 'bsxfun');

ใน code ข้างบน เราแสดงวิธีคูณให้ดู 3 แบบค่ะ
C1 คือ แบบ A*B ใน for loop ธรรมดา
C2 คือ การคูณตามนิยามการคูณเมตริกส์ในรูปแบบ vectorized
C3 คือ การคูณโดยใช้ bsxfun

จะเห็นได้ว่าการคูณแบบ bsxfun นั้นเร็วกว่าการคูณธรรมดาไปหลายขุมเลยค่ะ แต่ค่ะแต่ สมมติว่าเราเปลี่ยนขนาด matrix n จาก 3 เป็น 50 ล่ะ

พอขนาดเมตริกซ์ที่ต้องการคูณใหญ่ขึ้น การคูณแบบ vectorized กลับใช้เวลามากสุดและการคูณธรรมดาใช้เวลาน้อยสุด

ดังนั้นโพสนี้จึงขอสรุปว่า แม้ว่าการคูณเมตริกซ์สามารถทำได้หลายวิธี แต่ขนาดของเมตริกซ์ก็มีผลกับการเลือกใช้วิธีในการคูณเมตริกซ์นั้น ดังนั้นเลือกวิธีให้เหมาะกับงานค่ะ

ความเร็วในการคูณ 2D matrix กับ 3D array แบบ slice by slice

โพสนี้จะขอเปรียบเทียบความเร็วในการคูณ element-wise 3D matrix แบบ slice by slice อธิบายเป็นภาษาไทยลำบาก วาดรูปให้ดูเลยดีกว่า

สมมติเรามี 2D matrix A ที่อยากจะเอาไปคูณแบบ element by element กับ 2D matrix B ที่มีจำนวน r ชั้น เราสามารถเขียนโค้ดได้หลายแบบ แต่ในที่นี้จะเปรียบเทียบแค่ 2 แบบ คือแบบ vectorization กับ แบบใช้คำสั่ง bsxfun เราจะข้ามการคูณแบบใช้ for loop ไป เพราะรู้คำตอบอยู่แล้วว่าช้ากว่าทั้ง2แบบนี้มาก แต่เนื่องจาก MATLAB เค้าเป็นจ้าวแห่ง matrix operation (ชื่อเค้าก็บอกอยู่ MATrix LABoratory) และ bsxfun เป็น operation ที่ low level มากๆ เลยน่าสนใจว่าใครจะไวกว่ากัน 

โค้ดที่ใช้ในการเปรียบเทียบเป็นดังนี้

close all; clear; clc;

% 2D matrix times each slice of 3D matrix
N = [10,50,100,200,500,700];
t1 = zeros(1,length(N));
t2 = t1;
for ii = 1:length(N)
    n = N(ii);
    A  = rand(n,n);
    B  = rand(n,n,n);
    tic;
    C = bsxfun(@times, A, B);
    t1(ii) = toc;

    tic;
    D = ones(1,1,n);
    A = A.*D;
    C = A.*B;
    t2(ii) = toc;
end
h = plot(N,t1,'-*b',N,t2,'-*r');axis tight;
h(1).LineWidth = 2;
h(2).LineWidth = 2;
set(gca,'FontSize',14);
legend('bsxfun','vectorization');
xlabel('n'); ylabel('time(s)');

ผลสรุปออกมา คืออออออ...... bsxfun ชนะค่าาาาา ยิ่งmatrixมีขนาดใหญ่มากยิ่งทิ้งห่าง เพราะงั้นหากใครต้องการความไวในการเขียนโค้ด ดิชั้นแนะนำให้หันมาใช้ bsxfun จ้าา

imagesc vs pcolor vs surf

โพสนี้เราจะมาเล่าเกี่ยวกับคำสั่งที่ใช้พล๊อต 2D matrix กันค่ะ

ก่อนอื่นเราจะสร้าง 2D matrix array ขึ้นมาก่อนโดยใช้ตัวอย่างจากโพสก่อนๆก็แล้วกัน

%% imagesc vs pcolor vs surf
npts = 51;
x = linspace(-1,1,npts);
y = x;
[Y,X] = meshgrid(y,x);
[phi,rho] = cart2pol(X,Y);

เอาล่ะ สมมติว่าเราจะพล๊อต rho ซึ่งเป็น 2D array โดยปกติเราจะใช้คำสั่ง imagesc ค่ะ

%% imagesc
figure(1)
imagesc(rho); axis image;

ที่ใช้คำสั่งนี้เป็น default เพราะว่ากราฟจะแสดงค่าที่อยู่ในmatrixโดยตรง และสีก็จะscaleตามค่าของmatrixนั้น 

ข้อเสียของimagescคือ ถ้าเรามีจำนวน data ไม่มาก รูปที่ได้จะไม่smooth อย่างที่เห็นรูปข้างบนเราสามารถเห็นpixelของแต่ละช่องได้

ดังนั้นเราจึงใช้อีกคำสั่งคือ pcolor จะได้รูปแบบนี้

%% pcolor
figure(2)
pcolor(rho);axis image;

ถ้าใช้คำสั่งนี้โดยไม่เขียนอะไรเพิ่มเติม เราจะเห็นช่องตารางของแต่ละ pixel แต่ถ้าเราเพิ่มคำสั่งเข้าไปดังนี้

shading interp;

รูปที่ออกมาจะ smooth มากกว่าที่ใช้ imagesc

คำสั่งสุดท้ายที่จะพูดถึงคือ surf ซึ่งทำงานคล้ายกับ pcolor แต่ใช้ในการวาดรูป 3 มิติ

%% surf
figure(3)
surf(x,y,rho); axis image;

เช่นเดียวกับ pcolor ถ้า plot surf โดยไม่ได้ใส่คำสั่งเพิ่มเติม รูปที่ออกมาจะเห็นเส้นตารางคั่นแต่ละ pixel

ถ้าใส่

shading interp;

เส้นที่คั่นก็จะหายไป และสีก็จะ smooth ขึ้น

แต่นอกจากจะทำให้smoothได้แล้ว สำหรับ 3D plot เรายังใส่ lighting และ เลือกมุมที่จะดูได้ด้วย

camlight; lighting gouraud;
lightangle(10,90)
view(45,35);

รูปก็จะออกมาดูดียิ่งขึ้นไปอีก

ทำกราฟให้สวยขึ้นเพื่อใช้ในรายงาน

หลายครั้งที่เราพล๊อตกราฟแล้วตัวอักษรเล็ก เล้นกราฟเล็ก ดูยุ่งเหยิง แล้วถ้าจะเอาไปใส่ในรายงานก็ต้องเขียนกำกับแกนอีก วันนี้เลยมาแสดงตัวอย่างในการ edit graph ให้สวยขึ้นในการใช้ในรายงาน ถือว่าเป็น code ที่ใช้บ่อยมากเหมือนกัน

1. สร้างกราฟตัวอย่างก่อน วันนี้จะพล๊อต sinc function ก็แล้วกัน

%% graphic for report
x = linspace(-10,10,500);
y = sinc(x);

ถ้าใช้คำสั่ง plot(x,y) หน้าตาจะออกมาแบบนี้

เส้นก็บาง ตัวอักษรก็เล็ก แทบไม่เห็นอะไรเลยใช้มั้ย

มาเปลี่ยนหน้าตากัน

1.  เปลี่ยนสีเส้น ทำเส้นให้หนา และจัดหน้าต่างให้พอดีกับกราฟ

plot(x,y,'-b', 'LineWidth', 2);
axis tight;

2. ต่อมาก็ทำตัวหนังสือให้ใหญ่ขึ้น ทำกรอปให้หนาขึ้น เขียนกำกับแกน x,y และใส่ชื่อกราฟ ในตัวอย่างนี้จะใส่ชื่อกราฟแบบ 2 บรรทัดให้ดู โดยการใส่ string เข้าไปใน list ถ้าต้องการเขียนแค่บรรทัดเดียว ก็ใส่stringธรรมดา เช่น title('sinc function');

%% Set graphic view
set(gca, 'FontSize', 14, 'LineWidth',2);
xlabel('x-axis');
ylabel('y-axis');
title({'sinc', 'function'});

3. สุดท้าย ถ้าต้องการเปลี่ยน tick label ของแกน x ก็สามารถทำได้ดังนี้

%% set tick markings
xm = [-3*pi:pi:3*pi];
xt = {'-3\pi', '-2\pi','-\pi', '0' '\pi' '2\pi' '3\pi'};
set(gca,'XTick', xm, 'XTickLabel',xt)

xm คือ ตำแหน่งที่ต้องการใส่label
xt คือ ชื่อที่ต้องการให้แสดง

เสร็จแล้วค่ะ ... ทีนี้ก็ได้รูปที่ดูดีขึ้นไปแปะในรายงานแล้ว

วาดวงกลม จาก 1D vector

อีก 1 รูปที่ใช้บ่อย สมมติว่าเรามี vector อยู่ 1 อัน เราจะวาดวงกลมจาก vector นี้ได้อย่างไร

1. เหมือนเดิมค่ะ สร้าง index array ขึ้นมาก่อน

npts = 201;
x = linspace(-1,1,npts);
y = x;
[Y,X] = meshgrid(y,x);

2. เปลี่ยน coordinate x,y ให้เป็น polar coordinate โดย rho จะเป็นรัศมี และ phi คือมุม

[phi,rho] = cart2pol(X,Y);

phi จะมีค่าจาก -pi ถึง pi

3. สร้าง 1D vector ที่ต้องการจะทำเป็นวงกลมขึ้นมา ในที่นี้จะใช้ rand ฟังชั่นในการสร้าง โดยจำนวนจุดจะเท่ากับครึ่งหนึ่งของขนาด matrix

z = rand(1,ceil(npts/2));

4. หลังจากนั้นก็วาดรูปวงกลมด้วยคำสั่ง interp1 ตามข้างล่างนี้เลยค่ะ

Z = interp1(x(ceil(npts/2):end),z,rho);

ปล. คำสั่ง interp นี่ใช้ในการวาดรูปบ่อยเลยค่ะ มันดีจริงๆ

เพิ่มเติมอีกนิด ในส่วนนอกวงกลมจะมีค่าเป็น NaN (not a number) ซึ่งหลายๆครั้งใช้ในการคำนวณต่อไม่ได้ ดังนั้นก่อนจะนำรูปที่วาดไปเข้าสมการอื่น อาจจะเปลี่ยนNaN ให้เป็น 0 หรือเลขอื่นก่อนโดยใช้คำสั่งดังนี้

Z(isnan(Z)) = 0;

วาดวงกลมใน 2D matrix

นี่เป็นชุดคำสั่งที่ใช้บ่อยมากในงานผู้เขียน แต่กลับมีตัวอย่างให้เห็นน้อยมากในinternet ส่วนใหญ่จะวาดกันแบบเป็นเส้น แต่ในงานของผู้เขียนต้องการวาดวงกลมทึบในarrayเป็นส่วนใหญ่ วิธีการวาดก็ง่ายแสนง่าย

1. ตามเสต็ปค่ะ สร้าง index array ขึ้นมาก่อน ซึ่ง index array นี้ก็เปรียบง่ายๆคือการสร้างแกน x และ แกน y โปรดสังเกตุว่าในคำสั่ง meshgrid แกน y จะ อยู่หน้า แกน x เสมอ สลับกันไม่ได้นะคะ ในที่นี้เราจะให้แกน x,y มีขนาด 200 จุด

npts = 200;
x = linspace(-1,1,npts);
y = x;
[Y,X] = meshgrid(y,x);

2. สร้างวงกลม โดยในที่นี้กำหนดให้วงกลมมีรัศมีเท่ากับ 1

D = X.^2 + Y.^2 <= 1;

เท่านี้ก็เรียบร้อย ใช้คำสั่ง imagesc(D) ดูรูปได้เลย

หรือ ถ้าอยากจะสร้างวงรีก็ทำได้ไม่ยาก ตามนี้

npts = 200;
x = linspace(-1,1,npts);
y = 2.*x;
[Y,X] = meshgrid(y,x);
D = X.^2 + Y.^2 <= 1;
figure;
imagesc(D); axis image

Rotate matrix at arbitrarily angle

วันนี้มีปัญหาเกี่ยวกับการหมุน 2D เมทริกซ์ งมอยู่ซักพักก็คิดขึ้นมาได้ว่าทำยังไง

ก่อนอื่นต้องพูดว่า ถ้าอยากหมุนเมทริกซ์ 90 องศา มันมีคำสั่งง่าย คือ rot90(A) หรือง่ายกว่านั้นคือใช้คำสั่ง matrix transpose A.' แต่ถ้าอยากจะหมุน 45 องศาล่ะ ต้องทำยังไง 

สมมติ matrix ที่ต้องการหมุน, A, มีขนาด n x m

1. ก่อนอื่น สร้าง index array ขึ้นมาก่อน

x = linspace(-1,1,n);
y = linspace(-1,1,m);
[X,Y]=meshgrid(x,y);

2. กำหนดตัวแปร theta เพื่อบอกว่าจะหมุนกี่องศา และหมุน index array ตามสูตรนี้ (อ่านเพิ่มเติมเกี่ยวกับ rotation matrix ใน wikipedia ได้)

theta = pi/4;
Xq = X.*cos(theta) - Y.*sin(theta);
Yq = X.*sin(theta) + Y.*cos(theta);

3. หมุนmatrix A ด้วยคำสั่งinterp2

Vq = interp2(X,Y,A,Xq,Yq);

เท่านี้ก็เรียบร้อยแล้วค่ะ อยากจะหมุนกี่องศาก็เปลี่ยนตัวแปร theta เอา

หรือวิธีการหา Xq, Yq ที่สั้นกว่าแบบข้างบนคือ

[th, r] = cart2pol(X, Y);
[Xq, Yq] = pol2cart(th+pi/4, r);