Studio LGS™

안녕하세요 LGS입니다.
이번 시간은 저번 시간에 이어 FFT 계속입니다. MATLAB에서 직접 FFT 해보죠~
Fourier Analysis에 관해서 자세하게 알고 싶은 분은 아래 link를 참조하세요.
두 문서를 읽어 보시길 적극 추천합니다. 아주 쉽게 잘 설명해 놓았습니다.
Fourier Analysis 1
Fourier Analysis 2

다시 한번 말하지만 위의 사각파의 경우 수많은 Cosine의 합으로 표현 될 수 있습니다. FFT를 통해 얻고자 하는 우리의 궁극적인 목적은 어떤 주파수(fn)가  포함되어 있고 그 각각의 주파수의 크기(또는 기여도) An를 구하는 겁니다.

우선 간단한 Sine波 하나로 시작합시다.
Script를 하나 만듭니다. 아래 코드를 M-File Editor 창에 Copy & Paste !!
(코드는 blinkdagger.com를 참조하였습니다.)

fo = 4;   %frequency of the sine wave
Fs = 100; %sampling rate
Ts = 1/Fs; %sampling time interval
t = 0:Ts:1-Ts; %sampling period
n = length(t); %number of samples
y = 2*sin(2*pi*fo*t); %the sine curve
 
%plot the cosine curve in the time domain
sinePlot = figure;
plot(t,y)
xlabel('time (seconds)')
ylabel('y(t)')
title('Sample Sine Wave')
grid

Sampling Rate(Fs)란 1초에 몇번이나 데이터를 측정하겠나는 말입니다.
위의 코드 대로라면 1초에 100번 측정하는 군요. 즉 0.01초 마다 데이터를 측정한다는 말이죠(Ts)~
측정된 데이터 하나를 샘플(sample)이라고 합니다. 나머지는 문제없습니다.
F5를 살며시 눌러 줍시다. Sinewave.m 으로 저장하시구요~

FFT를 했을 때 우리가 기대하는 값은 4Hz (Hz = /초, 즉 1초에 4번 진동)와 크기 2입니다.
위의 Code를 보세요. 우리가 그린 sin 그래프는 fo(frequency)가 4이고 크기(Amplitude)가 2였습니다(y = 2*sin(2*pi*fo*t)). 그럼 FFT를 했을 때 4Hz 주파수에서만 Peak가 뜨고 나머지 주파수는 다 '0'이나와야 됩니다. 즉 위의 Sine은 4Hz 이외에는 어떤 주파수 성분도 담고 있지 않습니다.
자 그럼 FFT 해봅시다. MATLAB에서 FFT를 수행하는 함수는 fft 입니다.

>>fft(a)'
ans =
    .
    .
  -0.0000 +0.0000i

어라!! 주파수와 크기 정보를 얻으려고 했는데 이건 뭐여!?? 당췌~~ 냠냠...
그러나 걱정할 필요 없습니다.저만 잘 따라오세요~ ^^ 우선 아래 그림을 봐주세요
   

fft 함수를 사용했더니 복소수(Complex number)로 값을 돌려주는군요. 복소수는 크기와 위상(각도)정보를 담고 있습니다. 즉 fft 함수가 각 주파수의 Acos(2πf+θ)의 A와 θ의 정보를 복소수 형태로 돌려준것입니다.
지금 우리가 관심 있는건 바로바로바로 크기 A입니다. 어떤 주파수가 얼마나 들어있나에 관심있으므로 phase(위상)은 잠시 접어둡시다. MATLAB에서 복소수의 크기(Megnitude)를 돌려주는 함수는 abs 입니다.

그럼 크기 A를 그래프로 그려봅시다.
새로 Editor 창을 띄워서 아래 Code를 Copy & Paste!!


%plot the frequency spectrum using the MATLAB fft command
matlabFFT = figure;  %create a new figure
YfreqDomain = fft(y); %take the fft of our sin wave, y(t)
 
stem(abs(YfreqDomain));  %use abs command to get the magnitude
%similary, we would use angle command to get the phase plot!
%we'll discuss phase in another post though!
 
xlabel('Sample Number')
ylabel('Amplitude')
title('Using the Matlab fft command')
grid on
axis([0 100 0 120]) %resizing axis size

stem은 그래프의 한 종류입니다. Graph Catalog에서 확인해 보세요
F5 , FreqDomain.m 저장하시구요


오홋!! 뭔가 나온것 같기는한데... 그런데 자세히 보니 뭔가 좀 맞지 않군요 -_-;;
1. x축(Frequency)가 맞지 않습니다. 우리가 기대하는 값은 4Hz(Frequency)였는데 말이죠..
2. y축(Amplitude)도 맞질 않군요 -_-;; 흠.. 
3. 그리고 왜 Peak가 2개나 뜰까요?..
왜 위와 같은 결과가 나왔을 까요? 밑으로 내려가기 전에 한번 생각해보시기 바랍니다.

첫번째로 x축이 맞지 않는 이유는 우리가 FFT한 y에 시간 정보가 담겨져 있지 않기 때문입니다.
>>plot(y)해보세요. x축이 맞지 않는 이유와 같습니다. 그러므로 우리가 수동으로 x축을 수정해 주어야 합니다.
두번째로 y축을 Nomalize(%로 표현) 해주기 위해 sample수로 나주어 주면 됩니다.
마지막으로 쉽게 설명하자면 Cosine은 y축에 대칭이죠? cos(-1) = cos(1). 그래서 대칭된 2개의 Peak를 돌려줍니다. 즉 위의 그래프는 그래프의 w정중앙을 중심으로 서로 대칭입니다. 그럼 각각의 피크는 +4Hz와 -4Hz 이겠군요. 대개의 경우 + 주파수 값을 사용합니다. 대칭이니까 반은 잘라서 버려버리고 +주파수에만 집중합시다.
그런데 반으로 잘라 버리면 문제가 발생합니다. 예를 들어 2 = cos(2π) + cos(-2π) (주파수는 1Hz)인데 대칭이니까 cos(-2π)성분을 버려 버리면 1Hz의 기여도가 1/2로 줄어 들어 버립니다. 그래서 Normalize한 y값에 2를 곱해주어야 합니다.

이런 내용을 마음에 두고 다시 그래프를 그려 보겠습니다.
positiveFFT라는 함수를 하나 만들어 놓으면 두고 두고 편합니다.

function [X,freq]=positiveFFT(x,Fs)
   N=length(x); %get the number of points
   k=0:N-1;     %create a vector from 0 to N-1
  T=N/Fs;      %get the frequency interval
   freq=k/T;    %create the frequency range
   X=fft(x)/N*2; % normalize the data
 
   %only want the first half of the FFT, since it is redundant
   cutOff = ceil(N/2);
 
   %take only the first half of the spectrum
   X = X(1:cutOff);
   freq = freq(1:cutOff);

다른 것은 문제없고 ceil 이라는 함수만 알면 되겠군요. ceil함수는 올림을 해주는 함수입니다. 예를 들어 ceil(4.2)은 올림을 해서 5가 됩니다.

Command Window에 아래 코드를 붙여 넣어 보세요~

[YfreqDomain,frequencyRange] = positiveFFT(y,Fs);
stem(frequencyRange,abs(YfreqDomain));
xlabel('Freq (Hz)')
ylabel('Amplitude')
title('Using the positiveFFT function')
grid on
axis([0,20,0,1.5])

제대로 나왔군요~
위의 그래프를 분석해보자면 주파수는 4Hz이고 크기는 2가되겠습니다.
2*sin(2*pi*4*t)와 정확히 일치힙니다.
MATLAB FFT 끝~~ ^^


실전 연습해보죠~
제가 좋아하는 Radiohead 노래를 가지고 놀아 보죵 ^^

Radiohead의 Creep 시작부분에서 베이스 기타 소리만 죽여 보겠습니다~
FFT를 통해 주파수 분석한 다음 베이스 기타를 제외한 주파수만 살리고 베이스 기타 주파수는 죽여 버립시다.
그리고 다시 Time domain으로 되돌려 주면 됩니다.Frequency Domain에서 Time Domain으로 되돌리려면 ifft 함수(Inverse FFT)를 사용하시면 됩니다 그런데 베이스 기타 주파수가 어디냐구요?? 아무래도 베이스니까 주파수가 낮겠죠 그래서 낮은 주파수를 죽여 버렸습니다.

위의 프로그램은 제가 MATLAB GUI로 만든것입니다.(x,y축은 맞추질 않았습니다.)
원본 음악은 이렇습니다. 베이스 일렉 소리를 잘 들어보세요.


FFT를 통해 베이스 기타 소리를 제거 했습니다.
신기하죠? ㅋㅋ


이어폰으로 들을 때는 베이스 일렉이 잘 들렸는데 스피커로 들으니 잘 구분이 안가내요. -_-;
이번에는 심벌 소리만 뽑아 봤습니다.



zero pading 및 좀더 많은 정보를 원하시면 아래 링크 참조하세요~~
www.blinkdagger.com/fourier

오늘은 여기까지 입니다.
이제 FFT가 나와도 두려울 것이 없습니다. ㅋㅋ
다음 시간에는 공대생의 모르면 최대의 적이자 알면 가장 큰 무기!!!
미분방정식(Differential Equation)을 다룰 예정입니다.
그럼 다음 시간에 뵙겠습니다. 총총...

개인적인 사정으로 미분방정식 강좌는 무기한 연기되었음을 알려드립니다.
가능한 빠른 시일내에 여건이 허락하는 데로 강좌를 계속 진행하도록 노력하겠습니다.
이상입니다. -LGS