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 มันจะดีเสมอไป