MATLAB에서 1차 저역통과필터를 구현해보자

앗... 무쟈게 바쁜(^^) 와중에도 블로그의 스킨을 변경했습니다. ㅎㅎㅎ. 그러나 이전의 글들이 너무 주먹구구식으로 관리가 되다보니... 이전 글들은 틈나는 데로 조금씩 새로운 스킨에 맞게 바꾸어야겠어요....ㅠㅠ. 그래도 반응형 심플한 스킨을 적용해서 좋네요. 좀더 써보고 괜찮으면 결재해야죠^^ (사용 후 지급하는 유료형 스킨이라는..ㅠㅠ 그래도 이쁘니까요^^)

오늘은 요근래 좀 이야기했던 저역통과필터 중에서도 가장 간단한 1차 저역통과필터입니다. 이 이야기는 참 오래된 이야기인데요. 아~~~주 예전에 1차 저역통과필터(LPF)를 C로 구현하는 이야기를 먼저 헀거든요.

[The Robot/Prog.Lang.] - 저역통과필터와 고역통과필터를 C로 구현

그리고 나서.. 좀 잠잠이 있다가... 그 후에 난데없이 액셀에서 구현하는 저역통과필터와

[The Robot/Prog.Lang.] - 엑셀에서 저역통과필터 (Low Pass Filter) 구현하기

Python에서 구현하는 것까지

[The Robot/Prog.Lang.] - Python - IPython에서 구현하는 저역통과필터 Low Pass Filter

이야기를 했죠.. 결국.. 초간단 버젼 1차 저역통과필터를 C와 엑셀, Python으로 이야기를 했는데요. 흠.. 알고봤더니 MATLAB에서 이야기를 안했네요...^^ 사실 MATLAB 유저들은 MATLAB이 제공하는 엄~~청난 필터 디자인 툴 때문에 잘 안쓸것 같지만 뭐 그래도 살짝꿍 이야기를 할려구요^^.

이렇게 생긴 1차 저역통과필터를 m-file로 간단히 구현해볼려구요^^. 

일단.. 시험신호를 만들어야겠네요^^

% create raw date
sampleTime = 0.0001;
endTime = 2;
time = 0:sampleTime:endTime;

testY1 = 5*sin(2*pi*1*time);
testY2 = sin(2*pi*3*time);
testY3 = sin(2*pi*4*time);
testY4 = sin(2*pi*5*time);
testY5 = sin(2*pi*10*time);
testY6 = sin(2*pi*90*time);
testY7 = sin(2*pi*100*time);
testY8 = sin(2*pi*110*time);
testY9 = sin(2*pi*900*time);
testY10 = sin(2*pi*1000*time);
testY = testY1 + testY2 + testY3 + testY4 + testY5 + ...
    testY6 + testY7 + testY8 + testY9 + testY10;

figure; plot(time, testY); grid on
title('Raw Data')
set(gcf, 'Color', [1,1,1]);

이렇게 testY신호를 만들었습니다. 일단 진폭이 5로 제일 큰 주 신호는 1Hz짜리이구요. 그 뒤를 이어서 3, 4, 5, 10, 90, 100, 110, 900, 1000Hz를 더해서 살짝 노이즈 스럽게 신호를 만들었습니다. 이 아이의 그림은

위와 같이 생겼습니다. 약간 확대해서 보면

이렇죠^^. 이게 제가 의도한 데로 주파수 성분이 잘 있는지를 확인해보기 위해 FFT를 수행해보죠^^

% FFT raw data
sampleFreq = 1/sampleTime;
lenData = length(time);
NFFT = 2^nextpow2(lenData);
Y = fft(testY, NFFT)/lenData;
f = sampleFreq/2*linspace(0,1,NFFT/2+1);
figure; plot(f, 2*abs(Y(1:NFFT/2+1)));grid on
title('Raw Data FFT')
set(gcf, 'Color', [1,1,1]);

이렇게 코드로 짜면 알 수 있습니다. 아~~ 둘다 마직막에 set(gcf, 'Color', [1,1,1])이라고 한것은 캡쳐뜰때 하얀색 바탕화면이 보기 좋아서요^^ 아무튼 이 코드를 실행하면

요런 결과를 얻을 수 있습니다.

일단 900, 1000Hz 관촬되구요^^

90, 100, 110Hz도 보이네요^^

그리고... 1, 3, 4, 5, 10Hz도 다 있습니다. 사실.. 뭐 이걸 확인하는게 중요한 건 아닌데.ㅠㅠ.....

아무튼.. 이 수식을 m-file로 표현한 것은

% simple LPF
tau = 0.01;
LPFResult = zeros(1, length(time));
for i = 2:length(time)
    LPFResult(i) = (tau*LPFResult(i-1) + sampleTime*testY(i))/(tau + sampleTime);
end
figure; plot(time, testY, time, LPFResult, 'r'); grid on
title('LPF result')
set(gcf, 'Color', [1,1,1]);

입니다. 사실은 for문 안에 있는 한 줄이 전부이지요^^. 이렇게 tau를 0.01로 세팅하고 이제 그 결과를 보겠습니다.

아.. 그 전에 위와 같이 bode 선도를 그려본겁니다. 뭐 당연하지만.. 1차 LPF네요^^. 아~~ 저 보드 선도에서 살짝 특이한 것이 있는데요. 바로 x축의 단위가 Hz라는 겁니다. 디폴트는 MATLAB의 경우 rad/sec이거든요. 그걸 Hz로 바꾸어서 보기 편하게 한거죠^^

% draw bode
den = 1;
num = [tau, 1];
sys = tf(den, num);
h = bodeplot(sys); grid on
setoptions(h, 'FreqUnits', 'Hz', 'PhaseVisible', 'off');
eval(['title(''tau = ' num2str(tau) ' '');'])
set(gcf, 'Color', [1,1,1]);

그 부분은 위에 나타나 있습니다. Hz로 단위를 변경한 것은 setoptions에서 바꾸어 준거죠^^ 아.. 그래서 필터의 결과~~를 FFT를 수행한 것을 먼저 보면...

일단 900, 1000Hz 성분은 안보이네요.

그리고... 90, 100, 110Hz 성분은 많이 줄어든 (1 -> 0.2) 것이 확인됩니다.

실제로도 고주파 부분이 많이 제거되고 나타나네요^^

여기서 이제 만족할 수준으로 tau를 변경하면서 선정하면 되는 겁니다.^^. 다음에는 FDATOOL을 한 번 다뤼보고 싶네요^^ 아 그리고 위 그림을 얻은 전체 MATLAB 코드는

% create raw date
sampleTime = 0.0001;
endTime = 2;
time = 0:sampleTime:endTime;

testY1 = 5*sin(2*pi*1*time);
testY2 = sin(2*pi*3*time);
testY3 = sin(2*pi*4*time);
testY4 = sin(2*pi*5*time);
testY5 = sin(2*pi*10*time);
testY6 = sin(2*pi*90*time);
testY7 = sin(2*pi*100*time);
testY8 = sin(2*pi*110*time);
testY9 = sin(2*pi*900*time);
testY10 = sin(2*pi*1000*time);
testY = testY1 + testY2 + testY3 + testY4 + testY5 + ...
    testY6 + testY7 + testY8 + testY9 + testY10;

figure; plot(time, testY); grid on
title('Raw Data')
set(gcf, 'Color', [1,1,1]);

% FFT raw data
sampleFreq = 1/sampleTime;
lenData = length(time);
NFFT = 2^nextpow2(lenData);
Y = fft(testY, NFFT)/lenData;
f = sampleFreq/2*linspace(0,1,NFFT/2+1);
figure; plot(f, 2*abs(Y(1:NFFT/2+1)));grid on
title('Raw Data FFT')
set(gcf, 'Color', [1,1,1]);

% simple LPF
tau = 0.01;
LPFResult = zeros(1, length(time));
for i = 2:length(time)
    LPFResult(i) = (tau*LPFResult(i-1) + sampleTime*testY(i))/(tau + sampleTime);
end
figure; plot(time, testY, time, LPFResult, 'r'); grid on
title('LPF result')
set(gcf, 'Color', [1,1,1]);

% draw bode
den = 1;
num = [tau, 1];
sys = tf(den, num);
h = bodeplot(sys); grid on
setoptions(h, 'FreqUnits', 'Hz', 'PhaseVisible', 'off');
eval(['title(''tau = ' num2str(tau) ' '');'])
set(gcf, 'Color', [1,1,1]);

% FFT LPF
sampleFreq = 1/sampleTime;
lenData = length(time);
NFFT = 2^nextpow2(lenData);
Y = fft(LPFResult, NFFT)/lenData;
f = sampleFreq/2*linspace(0,1,NFFT/2+1);
figure; plot(f, 2*abs(Y(1:NFFT/2+1)));grid on
eval(['title(''FFT result of 1st LPF where tau = ' num2str(tau) ' '');'])
set(gcf, 'Color', [1,1,1]);

입니다.^^