Fourier Transform เบื้องต้น

ที่จริงการทำ Fourier Transform ใน MATLAB นั้นง่ายมาก แค่ใช้คำสั่ง fft ก็จบ แต่ส่วนใหญ่เท่าที่ทำ FFT มา จะมีปัญหาตลอดกับการทำ x-, y-axis ทำทีไร ลืมทุกที เพราะงั้น ก็จะมาจดโน๊ตไว้ ณ ที่นี้กันลืมแล้วกัน

ยกตัวอย่างง่ายๆ อย่างแรก เช่นการทำ Fourier transform ของ sine function

1. สร้างฟังชั่นตัวอย่างขึ้นมาก่อน

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1000;             % Length of signal 
t = linspace(0,(L-1)*T,L); % Time vector
f = 5;                % sine wave frequency
S = sin(2*pi*f*t);
plot(t,S); axis tight;
xlabel('time (s)'); ylabel('amplitude');

อ่ะ แล้วลงทำ fft ดู

Y = fft(S);
plot(abs(Y),'-b','LineWidth',2); axis tight;

ทำ fft แล้วจะได้กราฟด้านล่าง ซึ่งบอกไม่ได้ทั้งขนาดและความถี่ ดังนั้นเราจึงต้องมาทำแกนอีกที

เริ่มจากการทำขนาดแกน y ให้ถูกต้องก่อน

yaxis = fft(S./L);
plot(abs(yaxis),'-b','LineWidth',2); axis tight;

เมื่อ normalized แกน y แล้วจะเห็นได้ว่าขนาด peak นั้นจะลดลงครึ่งหนึ่งเสมอ ตามกฎของ fft

หลังจากนั้นมาทำแกน x ต่อ

NFFT = L;
xaxis = Fs/2*linspace(0,1,NFFT/2+1);
plot(xaxis,abs(yaxis(1:L/2+1)),'-b','LineWidth',2); axis tight;

อันนี้ตัดช่วงครึ่งหลังของ fft ออกไปนะ เพราะมันคือ mirror ของช่วงครึ่งแรก ก็จะได้กราฟตามนี้

ลองขยายกราฟดูว่าได้peakที่ตำแหน่ง 5 Hz จริงป่าว

เป๊ะเลยยยยยย

อ่อ NFFT คือจำนวนจุดที่ใช้ในการทำ fft ซึ่งบางทีอาจจะไม่เท่ากับ L ก็ได้ สมมติดังตัวอย่างต่อไปนี้

NFFT = 2^13;
yaxis = fft(S./L,NFFT);
xaxis = Fs/2*linspace(0,1,NFFT/2+1);
plot(xaxis,abs(yaxis(1:NFFT/2+1)),'-b','LineWidth',2); axis tight;

ตัวอย่างด้านบน เราทำ fft ไป 2^13 จุด เพราะอยากได้กราฟที่มีresolutionมากขึ้น ไม่เป็นแค่ peak เดียวอย่างกราฟข้างบนๆ ซึ่งในกรณีนี้มันคือ artifact แหละ เพราะfft ของ sine wave ที่ 1 frequency มันก็ต้องออกมา 1 frequency สิเนอะ แต่แค่ทำให้ดูเฉยๆว่ามันทำได้ แต่ไม่ได้หมายความว่าการเพิ่ม sampling rate มันจะดีเสมอไป

เริ่มการวิเคราะห์หุ้นด้วย 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 จ้าา