FDATOOL을 이용해서 Matlab에서 필터 설계하기

블로그 유입 키워드를 보면 matlab lowpass filter가 꽤 있다.

아마도 컴퓨터나 전기전자 공학 학생들인것 같고, 처음에 이런 필터를 어떻게 설계해야 하는지 찾다 찾다 여기 까지 왔으리라...



사실 이런 필터를 설계하는 방법은 여러가지가 있고, 직접 필터 설계하는 방법 및 공식을 코드로 작성할 수 도 있다.

그렇지만 매번 이렇게 코드를 작성하는 것은 너무 귀찮은 일이고 magnitude response, phase response, pole-zero plot등을 보기 위해서는 또 코드를 작성해야 할 것이다., 매틀랩에서는 다행히도 좋은 툴을 제공하고 있다.



커맨드 라인에서 fdatool을 입력해 보자. (단 signal processing toolbox가 설치되어 있어야 한다.)

왼쪽 아래 메뉴에 Response Type과 Design Method를 설정하는 부분이 있다.



Response Type은 어떠한 Response 특성을 같는 필터를 설계할 것인지 설정하는 부분이다.

Lowpass filter 뿐만 아니라, Highpass, Bandpass, Bandstop, Differentiator, Peaking, Notching, Multiband, Hilbert Transformer등 다양한 필터의 설계를 할 수 있다.


Design Method는 어떠한 방법으로 위의 Response Type을 가지는 필터를 설계할 것인지 설정하는 부분이다.

IIR, FIR을 설정할 수 있으며, 각 필터 타입 마다 여러가지 방법을 선택할 수 있다.



각 Design Method마다 특성이 어떤지는 구글링을 하던지 DSP 책을 참고하도록 하자.



일단 관심이 있는것이 Lowpass filter이므로 Lowpass filter를 선택하고, Design Method는 FIR에 window method 를 선택한다. 왜냐면 그냥 내 마음이다 ㅋㅋ (다른 이유가 있지만, 자세히 쓰기 귀찮으니 걍 넘어가자)

(그림엔 Equiripple이 선택되어 있지만 이것을 Window로 바꾸자)

여러가지 Window가 제공된다.웬만한 DSP 교재에 나와있는것은 다 들어있고, 맨 밑에 User Defined를 선택하면 사용자가 따로 만든 윈도우를 사용할 수 도 있다. 매틀랩에서는 이런 윈도우 함수를 만드는 툴(wintool)도 제공하고 있다.

나머지 정보를 입력하자.

Filter Order 부분은, 아주 쉽게 말하면 얼마나 긴 impulse response 를 가지는 필터를 설계할 것인지를 설정하는 것이다.

이것이 길면 길 수록, 필터의 실제 반응이 좋아지지만, 그만큼 많은 저장공간을 필요하고, 필터링에 의한 delay가 길어지게 된다.

또 너무 길어서 필터의 magnitude response가 rectangluar 형태에 가까워지면 Gibb's phenomenon에 의한 ripple이 신호에 크게 나타나므로 '적당히' 설정해야 한다.

또한, 실제로 DSP 보드에 필터를 설계해서 올린다면, 저장공간이나 필터에 의한 딜레이, quantization 또한 고려해야 한다.



밑에 비활성화 되어있는 Minimum order는 주어진 필터 스펙에서 최소 order를 갖도록 설정해주는 탭이다. 현재 Hamming Window방법은 제공하지 않는다. 다른 윈도우 방식에서는 적용되니까 본인들이 나중에 한번 해보자.


Frequency Specification부분은 Sampling frequency, Fs, 와 Cutoff Frequency, Fc, 를 설정하는 부분이다. 오른쪽에 Magnitude Specifications에 잘 나와있듯이 6dB 감소하는 부분을 의미한다.

Units라고 되어있는 팝업 메뉴는 주파수의 단위를 설정하는 부분이다. Normalized, Hz, kHz, MHz등을 선택할 수 있다.

이 메뉴들 - Frequency Specification, Magnitude Specification - 은 필터 설계 방법을 어떤 것을 선택하느냐에 따라서 메뉴가 생기기도 하고 사라지기 도 한다. 그때는 위의 Filter Specifications라고 되어있는 그림을 보면 무엇을 의미하는지 잘 알 수 있을것이다. 나중에 본인이 한번 해보길 바란다. Passband ripple, stop band gain등을 설정할 수 있다.

고민하지 말고 그냥 sampling frequncy는 1000Hz로 하고, cutoff frequency는 250Hz로 해보자.


그리고 마지막으로 밑에 있는 Design Filter라는 버튼을 누르면 좀 시간이 걸린 후 필터가 만들어지게 된다!!

자동으로 아래와 같이 Magnitude Response가 그려지게 된다.

그래프 위의 버튼들을 차례로 눌러보면 Magnitude Response, Phase Response, group delay,impulse response, pole-zero plot, coefficient, filter information등 필터 설계에 필요한 거의 모든 정보들을 볼 수 있다.

필터의 impulse response

필터의 pole-zero plot

필터의 coefficient

필터의 정보, 어떠한 타입으로 설계 되었는지, stability, 필요한 연산 수 등이 표시된다.

이제 설계된 필터를 사용하는 방법을 알아 보자. fdatool은 설계된 필터를 여러가지 포맷으로 export할 수 있다. 그냥 matlab의 workspace로 보낼 수도 있고, simulink object로 보낼 수도 있고, xilinx coefficient로 보낼 수 도 있다.

일단 그냥 matlab의 workspace로 보내보자.

File - Export 메뉴를 선택하자.

여러가지 메뉴가 또 나오는데, Export To는 어떠한 대상으로 필터를 출력할 것인지 선택하는 부분이다. Workspace로 할 것이므로 그냥 Workspace로 선택하자.

Export As는 어떠한 변수 형태로 출력할 것인지 선택하는 것이다. Coefficient만 할것인지 나머지 정보들까지 같이 Object로 출력할 것인지 선택한다. 여기 까지 쓰다보니 귀찮으니까 그냥 Coefficients로 선택하자.


Variable Names는 출력될 변수의 이름이다. 기본이 Num이다. (왜냐면 FIR 필터는 numerator만 있으니까 ㅋㅋ)


Export버튼을 누르자~

Workspace를 보면 Num이라는 1x21 크기의 행렬이 만들어져 있다. 이 안에 아까 보았던 impulse response가 들어있다.

이제 이 필터로 필터링을 해보자. 50Hz sine wave를 만들고, 여기에 AWGN (Additive White Gaussian Noise)를 넣고, 그것을 다시 필터링해서 노이즈를 제거할 것이다.


아래 처럼 명령어를 쳐서 50Hz sine wave를 만들어 보자.

만들어진 50Hz sine wave

15dB의 SNR을 가지도록 AWGN을 넣어주자. 다음과 같이 입력하자.

>> y2 = awgn(y, 15, 'measured');

아래는 그 결과이다.

이제 필터링을 하자.

필터링을 하는데는 여기서도 여러가지 방법이 있으나, 매틀랩의 filter라는 함수를 사용해 보자. filter 함수의 사용법은 매틀랩 help를 참고하자.

필터링 결과는 다음과 같다.

White Gaussian Noise는 모든 주파수 대역에 분포해 있으므로 우리가 만든 필터로는 완전히 노이즈를 제거할 수 없으나, 그래도 필터링 하기 전 보다는 훨씬 나아진 것을 알 수 있다.


그런데 왼쪽을 보면 그냥 평탄하게 신호가 없는 부분이 있는 것을 알 수 있다. 필터 딜레이에 의한 것이다.

FIR 필터는 이러한 딜레이를 없앨 수 있으며, 매트랩에서는 이것도 따로 함수를 제공한다. filtfilt함수 이다.

(Zero-Phase Filtering이라고 한다.)

filtfilt함수는 FIR 필터링에서 딜레이를 보정해준 후의 결과를 보여준다. 다음과 같이 입력하자.

아래는 filtfilt함수를 사용한 결과이다.

필터링이 잘 되었을 뿐 아니라, 딜레이까지 사라진 것을 볼 수 있다.





아마 짐작하고 있겠지만, 매틀랩에서는 fdatool을 이용한 설계를 그대로 매틀랩 코드로 출력할 수 있다.

File - Generate MATLAB Code 메뉴를 누르면 된다. haha.m으로 저장하였다.

function Hd = haha
%HAHA Returns a discrete-time filter object.

%
% MATLAB Code
% Generated by MATLAB(R) 7.13 and the DSP System Toolbox 8.1.
%
% Generated on: 26-Oct-2011 20:53:50
%

% Equiripple FIR Lowpass filter designed using the FIRPM function.

% All frequency values are in Hz.
Fs = 48000;  % Sampling Frequency

Fpass = 9600;            % Passband Frequency
Fstop = 12000;           % Stopband Frequency
Dpass = 0.057501127785;  % Passband Ripple
Dstop = 0.0001;          % Stopband Attenuation
dens  = 16;              % Density Factor

% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fpass, Fstop]/(Fs/2), [1 0], [Dpass, Dstop]);

% Calculate the coefficients using the FIRPM function.
b  = firpm(N, Fo, Ao, W, {dens});
Hd = dfilt.dffir(b);
% Set the arithmetic property.
set(Hd, 'Arithmetic', 'fixed', ...
    'CoeffWordLength', 16, ...
    'CoeffAutoScale', true, ...
    'Signed',         true, ...
    'InputWordLength', 16, ...
    'inputFracLength', 15, ...
    'FilterInternals',  'FullPrecision');
denormalize(Hd);


% [EOF]

이것은 다시말하면, 굳이 fdatool이라는 GUI를 통하지 않고서도, 위의 명령어를 이용해서 필터를 설계할 수 있다는 이야기이다. (사실 저 명령어 뿐만 아니라 다양한 필터 설계 함수가 존재한다. help를 뒤져보면 엄청 많다.)



공학인들, 쫄지말고 화이팅 하자.