간단히 MATLAB을 이용하여 체비세프 ChebyShev 저역통과 필터 구현해보기

사실... 저역통과필터는 뭐 원체 많이들 사용하고 있는거라 어려울게 없습니다만.... 그걸 저는 또 구지~~~ C Code로 구현하는 것에 대해 한 번 다룬적[바로가기]이 있었죠... 그리고는 그걸 다시 Python으로 구현하는 법을 이야기[바로가기]를 했구요. 또.. MATLAB으로 구현하는 것도 각각 FFT까지 수행해 가면서 다루었지요.[바로가기] 심지어는 엑셀에서 구현하는 방법마저도... 또 다루었습니다^^[바로가기] 이제 끝나나 하셨겠지만~~~^^ 이번에는 ChebyShev 체비세프 저역통과필터를 한 번 이야기해볼려고 합니다. 이 글을 쓰는 카테고리가 Program Language로... MATLAB에서 그냥 간편히 빠르게 어떻게 구현할 것인가가 목적입니다. 그러니 뭐 유도과정 원리... 등등의 복잡한 이야기는 [바로가기]로 패스~~~하고 시작하죠^^ 

일단 체비세프 필터의 주파수 특성은

출처 : 위키백과 Chebyshev filter 항목의 그림

와 같습니다. 통과를 시키는 저역부분에 리플(ripple)이 존재한다는 단점이 있지만...

출처 : 위키백과의 Chebyshev필터

Butterworth에 비해 주파수를 자르는(cut-off) 기울기가 아주 가파르게 잘 자르는 특성을 알 수 있습니다. 오늘 다루는 것이 체비세프 1형이고.. 2형의 경우는 통과대역의 ripple은 없지만.. 차단된 대역에 ripple이 있다는 단점이 있습니다.

체비세프 필터를 MATLAB에서 설계하는 간편한 절차는 위와 같습니다. passband와 stopband를 정하고... ripple의 정도를 설정한 후~~ cheb1ord를 실행하면 필터의 차수와 가중치를 계산해 줍니다. 그것과 통과대역의 ripple정도를 설정한 후에 cheby1을 실행하면 필터의 계수를 정해 줍니다. 이를 filter라는 함수를 통해 원 신호에 필터를 적용한 결과를 알 수 있게 해 줍니다. filter라는 함수를 나중에 C든... 뭐든 코드로 구현하면 되는데요... 그 원리는

출처 : MATLAB filter 명령의 설명 문서

입니다. 이산 시간에서의 개념이므로 차수에 맞게 이전 값을 저장해서 계산해 주면 됩니다. 일단 오늘은 MATLAB~~~이니까.. 예제 전체를 설명해야죠~~^^. 더불어 [바로가기]에서 했던 1차 LPF를 가지고 비교해볼까 합니다.^^

% create rawData
sampleTime = 0.0001;
endTime = 2;
time = 0:sampleTime:endTime;

freq = [1 2 3 4 10 70 80 90 100 200 300 400 500 1000 1100 1200 1300 2000];

testYtmp = zeros(length(freq), length(time));

for i=1:length(freq)
    eval(['testYtmp(' num2str(i) ', :) = sin(2*pi*freq(i)*time); '])
end
testYtmp(1,:) = testYtmp(1,:)*5;

refSignal = testYtmp(1, :) + testYtmp(2, :) + testYtmp(3, :) + testYtmp(4, :) + testYtmp(5, :);
testY = sum(testYtmp, 1);

이 부분은 시험 신호를 만드는 부분입니다. 먼저 1, 2 Hz 대역의 신호 몇 개와 10 70 등등의 신호들... 수백 Hz 신호들을 조합해서 test 신호를 만들었습니다. 제 목적은 10Hz까지는 확실히 복원하고 나머지는 버리는 필터를 설계하고 싶은거구요. 그래서 refSignal이라는 것에서 10Hz까지를 따로 두고 나중에 비교할려고 합니다.

%% simple 1st LPF
tau = 0.018;
resultLPF = zeros(1, length(time));
for i = 2:length(time)
    resultLPF(i) = (tau*resultLPF(i-1) + sampleTime*testY(i))/(tau + sampleTime);
end

이 부분은 이전에도 이야기한 1차 LPF를 적용한 결과이구요.... 저~~ 필터 계수 tau는 오늘 세팅할 체비세프 필터와 응답 시간을 최대한 맞출려고 그냥 막 잡은 겁니다. 비교를 위해서이죠~~~^^

%% ChebyShev LPF
passband = 0.01;
stopband = 0.029;
passrip = -20*log10(0.99);
stopatten = -20*log10(0.0001);

[Nc1, Wnc1] = cheb1ord(passband, stopband, passrip, stopatten);
[Bc1, Ac1] = cheby1(Nc1, passrip, Wnc1);

resultCheby = filter(Bc1, Ac1, testY);

오늘의 핵심인 체비세프 필터입니다. 위에서 설명했습니다만~~~^^

%% Draw the result
figure; plot(time, testY, 'LineWidth', 0.1)
grid on
hold on
plot(time, refSignal, 'c', 'LineWidth', 1.5)
plot(time, resultLPF, 'k', 'LineWidth', 1.5)
plot(time, resultCheby, 'r', 'LineWidth', 1.5)
hold off
set(gcf, 'Color', [1,1,1])
legend('rawData', 'refSignal', '1stLPF', 'ChebyShev')
eval(['title(''1st LPF coef = ' num2str(tau) ', ChebyShev''''s order = ' num2str(Nc1) ' '');'])

이제 원신호(testY)와 1차 LPF(resultLPF)와 체비세프 필터(resultCheby)를 비교하기 위해 그래프로 그려보죠~~^^

위 결과와 같습니다. 파란색은 노이즈를 잔뜩 실어준 원 데이터이구요. 하늘색은 딱 저 모양으로 복원되었으면 하는 바램(^^)이구요. 까망색은 1차 LPF, 빨간색이 체비세프 필터의 결과입니다. 두 필터는 성능 비교를 위해 지연된 시간을 대충 눈으로(^^) 일치시켜서 보는겁니다.

이 것만 봐도 체비세프가 잔 진동 없이.. 이말은 다른 주파수 영역없이 상대적으로 깔끔하게 잘 필터링 된 것이라고 볼 수 있는데요...

%% FFT
sampleFreq = 1/sampleTime;
lenData = length(time);
NFFT = 2^nextpow2(lenData);

Y1 = fft(testY, NFFT)/lenData;
Y2 = fft(resultLPF, NFFT)/lenData;
Y3 = fft(resultCheby, NFFT)/lenData;

f = sampleFreq/2*linspace(0,1,NFFT/2+1);

figure; plot(f, 2*abs(Y1(1:NFFT/2+1)), ...
    f, 2*abs(Y2(1:NFFT/2+1)), 'c', ...
    f, 2*abs(Y3(1:NFFT/2+1)), 'r', 'LineWidth', 1.5 )
grid on
set(gcf, 'Color', [1,1,1])
legend('rawData', '1stLPF', 'ChebyShev')

좀더 정확히 알기 위해 FFT를 수행해보죠. 위 코드는 뭐 이미 제가 이야기한 적이 있는 FFT코드이구요.. 

결과가 위와 같습니다. 원신호는 제가 의도한 것 처럼 주파수가 잘 버무려진듯 보이네요...

확실히 하늘색(1차 LPF)는 점점 주파수 성분이 깍여가는 것이 보이네요. 좀 더 보면

1차 LPF는 자르고 싶은 주파수 대역이 확실히 깔끔하게 자르지 못하고 남은 영역이 보이지만. 빨간색 체비세프 필터는 상대적으로 잘 잘려진것이 확인됩니다.^^

다시 원 신호를 더 자세히 보면... 이 의미를 확인할 수 있습니다.^^. 잘 잘라서 10Hz까지만 유지를 하고.. 나머지는 잘 잘렸네요~~~^^