Module 9. FFT Filter (Fast Fourier Transform) - (4)Gaussian 필터

FFT 필터 관련 첫글에 첨부해두었던 mcrFFT.bas를 엑셀에서 열어보시면 몇 개의 함수로 구성되어있습니다. 함수명을 보시면 대략적인 기능을 알 수 있을 것입니다.

N개의 실수로 구성된 데이터를 Fourier 변환하기 위해서는 FFTR1D 함수를 사용하면 됩니다.

인덱스가 (0~N-1)인 데이터셋을 입력하여 FFT를 실행하게 되면, (0~N-1)의 인덱스를 갖는 복소수 배열이 출력됩니다.

동일한 방법으로 (0~N-1)의 인덱스를 갖는 복소수를 FFTR1DInv 함수로 역변환하게 되면, 실수 배열을 출력해줍니다. 따라서, 전체적인 계산 순서를 정리하면 아래와 같습니다.

실수 배열(입력) → (FFTR1D) → 복소수 배열 → (Filter 적용) → (FFTR1DInv) → 실수 배열 (출력)

FFTR1D 함수와 FFTR1DInv 함수는 위의 모듈에 이미 만들어져 있으므로, 데이터 형식에 맞게 Filter를 어떻게 만들 것인지만 고민하면 됩니다.

일반적으로 많이 사용되는 Gaussian 필터의 경우, 아래와 같이 작성할 수 있습니다.

---------------------------------------------------------------

'Highpass 또는 Low Pass 필터인 경우

Public Function FFT_GaussianFiltered(iData() As Double, iCutOff As Double, Optional iHighpass As Boolean = False) As Double()
  Dim tFFT() As typeXYPair, tFFTInv() As Double, tNC As Double, tVal As Double, n As Long, i As Long

  '입력받는 데이터는 0~n까지 index를 가짐
  n = UBound(iData)

  '실수값을 갖는 1차원 배열을 FFT 실행하여 복소수 형식의 데이터를 출력받음
  FFTR1D iData, n + 1, tFFT

  'CutOff값을 기준으로 공통으로 계산이 반복되는 부분만 미리 계산하여 변수이 저장해둠 (단순히 계산 횟수를 줄이기 위한 방법)
  tNC = iCutOff / n

  'High Pass인 경우, 높은 주파수(i가 증가하는 경우)에서는 투과율이 1에 수렴하고, Cutoff에서는 0.5, 낮은 주파수에서는 0에 수렴하는 filter값을 가지며, low pass인 경우에는 반대가 됨
  Select Case iHighpass
    Case True
      For i = 0 To n
        tVal = 1 - 2 ^ (-(i * tNC) ^ 2)
        tFFT(i).X = tFFT(i).X * tVal
        tFFT(i).Y = tFFT(i).Y * tVal
      Next
    Case False
      For i = 0 To n
        tVal = 2 ^ (-(i * tNC) ^ 2)
        tFFT(i).X = tFFT(i).X * tVal
        tFFT(i).Y = tFFT(i).Y * tVal
      Next
  End Select

  '필터를 곱한 결과를 FFT 역변환을 실시합니다.

  FFTR1DInv tFFT, n + 1, tFFTInv
  FFT_GaussianFiltered = FFTR1DInv
  Erase tFFT, tFFTInv
End Function

---------------------------------------------------------------

'Band pass 필터인 경우, 낮은 주파수의 cutoff값과 높은 주파수의 cutoff값을 입력받은 후, 두 cutoff 사이의 영역만 투과시키는 filter가 됨
Public Function FFT_GaussianBandFiltered(iData() As Double, iCutOff1 As Double, iCutOff2 As Double) As Double()
  Dim tFFT() As typeXYPair, tFFTInv() As Double, tNS As Double, tNC As Double, tVal As Double, n As Long, i As Long
  n = UBound(iData)
  FFTR1D iData, n + 1, tFFT
  tNS = iCutOff1 / n
  tNC = iCutOff2 / n
  For i = 0 To n
    tVal = 2 ^ (-(i * tNS) ^ 2) * (1 - 2 ^ (-(i * tNC) ^ 2))
    tFFT(i).X = tFFT(i).X * tVal
    tFFT(i).Y = tFFT(i).Y * tVal
  Next
  FFTR1DInv tFFT, n + 1, tFFTInv
  FFT_GaussianBandFiltered = FFTR1DInv
  Erase tFFT, tFFTInv
End Function

---------------------------------------------------------------

위의 계산에서 Cut-Off라는 용어를 계속 사용하였습니다. 여러가지 측정이나 계산목적으로 FFT 변환을 사용해보신 분들은 cut-off라는 용어가 익숙하실 수 있습니다만, 일반인들이라면 cut-off, 주파수... 등의 개념이 익숙하지 않을 수 있기에 약간의 부연 설명을 드립니다.

시간 함수를 기준으로 설명을 드리면, 예를 들어, 1초 간격으로 수집한 데이터가 있다고 하겠습니다. 이때, 데이터는 0초에 측정된 데이터, 1초에 측정된 데이터, ..., n초에 측정된 데이터로 구성됩니다. 만약, 이 데이터를 FFT 변환하여 Cos, Sin 함수로 분해해주었다고 하겠습니다.

FFT 변환을 하게 되면, 변환결과는 전체 함수의 평균값(주파수가 0), 주파수가 1인 Cos, Sin 함수, 주파수가 2인 Cos, Sin 함수, ...., 주파수가 n인 Cos, Sin 함수로 나눠지게 되며, 각 함수들의 진폭을 데이터 배열에 할당해줍니다. 따라서, index가 낮은 데이터는 낮은 주파수, index가 높은 데이터는 높은 주파수의 진폭값이 됩니다. (여기에서 주파수가 1, 2, 3,... 이라고 하는 의미는 측정된 전체 구간을 주기로 하는 함수, 1/2 구간을 주기로 하는 함수, 1/3 구간을 주기로 하는 함수.... 를 의미합니다.)

만약, 주파수가 10 이상이 되는 Cos, Sin 함수는 제거(filter=0)하고, 낮은 주파수의 함수만 다시 결합하게 되면, 전체적으로 노이즈가 제거된 형태의 함수를 얻게 되는데, 이러한 필터를 low pass filter라고 하고, 반대로 high pass filter의 경우에는 높은 주파수의 함수만으로 조합되므로 노이즈만 남고 background가 제거되는 효과가 있습니다.

비록 주파수는 시간의 역수 단위, 주기는 시간 단위를 갖게 되는 것입니다. 주파수가 비록 익숙한 용어이기는 하나, 1초 단위로 측정된 데이터를 주파수 0.1Hz를 기준으로 필터링한다고 하면, 쉽게 그 의미가 와닿지 않습니다. 이때, high pass / low pass filter의 기준값이 되는 주파수에 해당하는 주기를 cut-off라고 부릅니다. 1초 단위로 측정된 데이터를 10초의 주기를 갖는 cos, sin 함수를 기준으로 장주기 데이터만 걸러 낸다거나, 단주기 데이터만 걸러낸다고 하면 조금은 상상하기 쉽지요.

다시말해, high pass/low pass 함수를 작성할 때, 입력값을 주파수, 혹은 파수번호로 입력하지 않고 cut-off를 입력으로 해두는 이유는 데이터의 X변수와 필터링하려는 기준인 cut-off가 동일한 단위를 갖는 값이기 때문에 직관적으로 이해하기 편하기 때문입니다.

위의 Gaussian 필터 함수에서 데이터를 입력한 후, cut-off를 10으로 지정한다면, 10개의 데이터를 주기로 갖는 cos, sin 함수를 기준으로 장주기 혹은 단주기 영역을 제거해주는 역할을 하게됩니다.