바로 프로젝트 적용 가능한 FIR Filter (low/high/band pass filter )를 c나 python으로 만들기

 band pass filter 코드가 필요해 인터넷을 검색을 해봤지만, 쓸 만한 코드와 자료가 별로 없어 그냥 만들어 쓰기로 했다.

필요한 공식은 위키피디아에 나와있다.

https://en.wikipedia.org/wiki/Finite_impulse_response

1. filter 크기

FIR filter를 만들 때 filter의 크기(tpas의 크기)를 결정할 때 도움될 수식이 있다.

  N=A㏈Fs / 22*Δf

A㏈ : 감쇄될 db
Fs : Sampling Frequency
Δf : frequency bandwidth

출처 : https://www.allaboutcircuits.com/technical-articles/design-of-fir-filters-design-octave-matlab/

def estimatefilterlen(fl,fh,fs,db):
    N = int(np.round(db * fs / (22 * (fh-fl)))-1)
    return N

2. Low Pass Filter

[ python code ]

def lowpassfilter(f,fs,N):
    taps = np.zeros(N)
    
    fc = f / fs
    omega = 2*np.pi*fc

    middle = int(N/2)
    for i in range(N):
        if i == middle:
            taps[i] = 2*fc
        else:
            taps[i] = np.sin(omega*(i-middle))/(np.pi*(i-middle))
        
    return taps

[ c code ]

#ifndef PI
#define PI 3.14159265358979323846
#endif

double *lowpassfilter(double f, double fs, int N)
{
    double *filter = NULL;
    double fc = 0;
    double omega = 0;
    int i = 0;
    int middle = 0;

    if (f <= 0 || fs <= 0 || N <= 0)
        return NULL;

    filter = (double*)malloc(N*sizeof(double));
    if (!filter)
        return NULL;

    memset(filter, 0, N*sizeof(double));

    middle = (int)(N / 2);
    fc = f / fs;
    omega = 2 * PI*fc;

    for (i = 0; i < N; i++)
    {
        if (i == middle)
        {
            filter[i] = 2 * fc;
        }
        else
        {
            filter[i] = sin(omega*(i - middle)) / (PI*(i - middle));
        }
    }
    return filter;
}

f = 2.0, fs = 100, N=908 인 low pass filter는 아래와 같은 모양과 주파수 응답 특성을 지닌다.

만약 위 1번의 필터 크기 결정 공식에 따라 필터의 크기를 달리 했을 경우의 예는 아래와 같다.

필터의 크기에 따라 주파수 응답 특성이 다르게 나타나는 것을 볼 수 있다. 필터의 길이가 길수록 주파수 응답 특성이 좋아지나 실시간 필터링의 경우 N/2 만큼의 딜레이가 발생하기 때문에 주의 깊게 사용해야 한다.

3. High Pass Filter

[ python code ]

def highpassfilter(f,fs,N):
    taps = np.zeros(N)
    
    fc = f / fs
    omega = 2*np.pi*fc
    
    middle = int(N/2)
    for i in range(N):
        if i == middle:
            taps[i] = 1 - 2*fc
        else:
            taps[i] = - np.sin(omega*(i-middle))/(np.pi*(i-middle))
                
    return taps

[ c code ]

double *highpassfilter(double f, double fs, int N)
{
    double *filter = NULL;
    double fc = 0;
    double omega = 0;
    int i = 0;
    int middle = 0;

   if (f <= 0 || fs <= 0 || N <= 0)
       return NULL;

    filter = (double*)malloc(N*sizeof(double));
    if (!filter)
        return NULL;

    memset(filter, 0, N*sizeof(double));

    middle = (int)(N / 2);
    fc = f / fs;
    omega = 2 * PI*fc;

    for (i = 0; i < N; i++)
    {
        if (i == middle)
        {
           filter[i] = 1 - 2 * fc;
        }
        else
        {
            filter[i] = - sin(omega*(i - middle)) / (PI*(i - middle));
        }
    }
    return filter;
}

f = 2.0, fs = 100, N=908 인 high pass filter는 아래와 같은 모양과 주파수 응답 특성을 지닌다.

4. Band Pass Filter

[ python code ]

def bandpassfilter(fl,fh,fs,N):
    taps = np.zeros(N)
    
    flc = fl/fs
    fhc = fh/fs
    omegal = 2*np.pi*flc
    omegah = 2*np.pi*fhc
    
    middle = int(N/2)
    for i in range(N):
        if i == middle:
            taps[i] = 2*fhc - 2*flc
        else:
            taps[i] = np.sin(omegah*(i-middle))/(np.pi*(i-middle)) - np.sin(omegal*(i-middle))/(np.pi*(i-middle))
    
    return taps

[ c code ]

double *bandpassfilter(double fl, double fh, double fs, int N)
{
    double *filter = NULL;
    double flc = 0;
    double omegal = 0;
    double fhc = 0;
    double omegah = 0;
    int i = 0;
    int middle = 0;

    if (fl <= 0 || fh <= 0 || fl >=fh || fs <= 0 || N <= 0)
        return NULL;

    filter = (double*)malloc(N*sizeof(double));
    if (!filter)
        return NULL;

    memset(filter, 0, N*sizeof(double));

    middle = (int)(N / 2);
    flc = fl / fs;
    omegal = 2 * PI*flc;
    fhc = fh / fs;
    omegah = 2 * PI*fhc;

    for (i = 0; i < N; i++)
    {
        if (i == middle)
        {
            filter[i] = 2*fhc - 2 * flc;
        }
        else
        {
            filter[i] = sin(omegah*(i - middle)) / (PI*(i - middle)) - sin(omegal*(i - middle)) / (PI*(i - middle));
        }
    }
    return filter;
}

fl = 1.0, fh = 2.0, fs = 100, N=908 인 band pass filter는 아래와 같은 모양과 주파수 응답 특성을 지닌다.

5. filtering

만들어진 FIR Filter를 사용하여 실제 샘플데이터를 필터링하는 코드는 아래와 같다. python코드가 간단해 python코드로 설명하겠다.

def filtering(samples, taps):
    N = taps.shape[0]
    out = np.zeros_like(samples)
    
    if len(samples) <= N:
        return out
        
    buffer = [samples[0] for _ in range(N)]
    buffer = np.hstack((buffer,samples))
    
    for i in range(len(samples)):
        out[i] = np.sum(buffer[i:i+N]*taps.T)
    
    return out[int(N/2):]

필터링은 convolution 연산을 사용한다. 필터의 길이에 해당하는 샘플과 필터의 합성곱의 결과들이 필터링 된 데이터가 된다. 연산을 위해 샘플 앞쪽에 필터의 길이만큼의 빈 데이터나 혹은 샘플의 끝부분 데이터 등을 넣어 추가해야 한다. 

필터링의 예제는 아래와 같다.

위와 같은 신호 데이터에 low pass filter를 적용하면 아래와 같은 결과를 얻을 수 있다.

6. windowing ***

FIR Filter에 윈도우를 적용하면 주파수 특성이 좋아지는 것을 볼 수 있다. 아래는 Band Pass Filter에 Hamming 윈도우를 적용해 비교한 것이다. 

    bpf = bandpassfilter(fl=fl,fh=fh, fs=fs,N=N)
    bpfw = bandpassfilter(fl=fl,fh=fh, fs=fs,N=N)*np.hamming(N)

주파수 특성을 보면 잘라져 나가는 주파수 영역이 적용하지 전에 비해 깨끗해지는 것을 볼 수 있다. 

실제 샘플에 적용해 보면 아래와 같은 차이를 볼 수 있다.

위와 같이 단순 Band Pass Filter를 적용할 경우 저주파 영역에 주파수 성분이 나타나는 것을 볼 수 있으나, 아래처럼 hamming윈도우를 filter에 적용하여 필터링 할 경우 저주파 영역의 주파수 성분이 나타나지 않는다. 

관련 글

FFT용 윈도우 함수 구현 예제 소스 코드