FFT를 이용한 noise 제거하기

nosie 제거 방법

신호를 측정시 신호에는 측정 대상이 되는 신호이외의 여러 noise들이 끼어 들어오게 된다. 

이는 원래의 신호를 분석하는데 noise들이 방해가 되므로 이들을 제거해야한다.

 

신호를 측정시 측정된 신호간에는 측정기기로 인한 term이 생겨 이산적인 데이터를 얻게 되므로

이를 푸리에 변환하기 위해서는 DFT를 이용해하지만 해당 방법은 연산에 드는 자원이 많이 든다.

 

DFT와 같이 이산적 데이터들을 푸리에 변환하면서 연산의 횟수를 줄이는 방법인 FFT를 사용한다.

데이터들을 FFT하여 시간 domain에서 진동수 domain으로 바꾸면 각 데이터점들이 어떠한 진동수들을 가지고 있는지 알 수 있다.

 

이때 측정하려는 신호는 일정한 진동수들의 조합으로 되어있으므로 각 데이터들점들은 해당 진동수들을 모두 포함하고 있을 것이다. 

임의의 noise들이 각 데이터 점에 껴들어가므로 원래의 신호가 가지고 있는 진동수와 같이 모든 데이터점에 포함되어 있지는 못한다.

따라서 nosie가 섞인 신호를 받아도 이를 진동수 영역으로 넘겨서 분석해보면 원래 신호가 갖고 있는 진동수가 다른 진동수들에 비해 더 클 것이다.

이를 이용하여 일정 크기에 미치지 못하는 진동수들은 모두 제거해버리고 남은 진동수만 가지고 푸리에 역변환을 하게 되면 측정하려는 원래 신호가 가지고 있던 진동수로 구성된 파형을 얻을 수 있다.

 

 

 

- 원래의 신호가 진동수 50과 120으로 구성되어 있으며 이 신호에 임의의 noise를 주었다.

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [16, 12]
plt.rcParams.update({'font.size': 18})

# 데이터간의 간격과 정의역 정의

dt = 0.001                                      # 각 데이터를 sampling 하는 간격
t = np.arange(0,1,dt)                           # 전체 데이터가 나타난 정의역 

# 잡음이 없을 때의 신호의 파형 

f = np.sin(2*np.pi*50*t) + np.sin(2*np.pi*120*t) # 진동수 50과 120을 갖는 파형이 합쳐져 만든 원래 신호
f_clean = f                                      # 원래신호를 clean 함수로 변수를 변경해두고 

# 원래 신호에 noise 추가

f = f + 2.5*np.random.randn(len(t))              # t개의 원래 데이터점 f(t)에 대해서 임의의 값을 2.5배 스케일링하여 더한다.

plt.plot(t,f,color='c',linewidth=1.5,label='Nosiy')
plt.plot(t,f_clean,color='k',linewidth=2,label='Clean')
plt.xlim(t[0],t[-1])
plt.legend()

 

검정색으로 나타난것이 원래 신호이며 빨간색으로 나타난것이 noise가 섞여 들어간 신호이다. 

우리는 검은색으로 나타난 신호를 원하지만 실제 측정을 통해 얻은 신호은 빨간색으로 나타난 신호이다.

따라서 FFT를 이용하여 빨간색 신호를 검은색 신호로 바꾸는 작업이 필요하다. 

 

 

n = len(t)                                 # FFT를 수행할 데이터의 점의 개수
fhat = np.fft.fft(f,n)                     # noise가 포함된 신호에 대해 FFT 수행

PSD = fhat * np.conj(fhat) / n             # Power spectrum (power per freq) 진동수 밀도
                                           # (a+bi)(a-bi)=a^2+b^2이므로 각 fhat에 conjugate를 곱한뒤 /n으로 정규화하면
                                           # 각 fhat의 크기를 얻을 수 있다. 
    
freq = (1/(dt*n)) * np.arange(n)           # (1/주기 = 진동수) 정의역을 시간에서 진동수로 바꾸어줌
L = np.arange(1,np.floor(n/2),dtype='int') # 진동수 정의역이 너무 크므로 반만 보기위해 설정

fig,axs = plt.subplots(2,1)

plt.sca(axs[0])
plt.plot(t,f,color='r',linewidth=1.5,label='Noisy')
plt.plot(t,f_clean,color='k',linewidth=2,label='Clean')
plt.xlim(t[0],t[-1])
plt.legend()

plt.sca(axs[1])
plt.plot(freq[L],PSD[L],color='c',linewidth=2,label='Nosiy')
plt.xlim(freq[L[0]],freq[L[-1]])
plt.legend()
plt.show

위의 빨간색으로 나타난 잡음이 섞인 시간 domain 신호를 FFT를 이용해 푸리에변환을 하여 진동수 domain으로 나타낸 것이다.

우리가 처음 측정하고자 했던 신호는 진동수 50과 120으로 구성되어 있는 신호였다.

위의 파란색 그래프에서도 50과 120부근에서의 진동수 밀도가 크게 나타나는 것을 알 수 있다.

측정된 데이터들이 원래 측정하려는 파의 진동수를 모두 가지고 있기 때문에 다른 noise들에 비해 더 높은 진동수 밀도를 갖는 것이다.

 

## PSD값을 이용하여 noise를 fittering

indices = PSD > 100       # PSD값이 100이상인 경우만 True 나머지 모두 False
PSDclean = PSD * indices  # PSD값에 indices를 곱하면 PSD값이 100이상인 진동수의 PSD값만 남는다.
fhat = indices * fhat     # noise_fhat에 indices를 곱하면 PSD값이 100이상인 진동수의 푸리에계수만 남는다. 
ffilt = np.fft.ifft(fhat) # 남은 푸리에계수만 푸리에 역변환을 하면 noise가 fittering된 파형을 얻게 된다. 

plt.plot(t,f_clean,color='k',linewidth=1.5,label='Clean')
plt.plot(t,ffilt,color='b',linewidth=2,label='Filtered')
plt.xlim(t[0],t[-1])
plt.legend()

검은색이 위에서 봤던 진동수 50과 120으로 구성된 원래의 파형이고 파란색이 PSD를 이용해 fittering된 푸리에계수를 역변환한 결과다.

그림에서 나타난것과 같이 거의 일치하는 것을 볼 수 있다. 따라서 FFT를 이용하면 noise를 제거하여 원래의 신호를 얻을 수 있게 된다.

 

 

Fittering하기전의 푸리에계수(각 진동수의 밀도,PSD)가 빨간색이고 Fittering이후가 파란색이다. 

Fittering을 통해 원래 파형을 구성하고 있는 진동수의 푸리에계수만 남기고 나머지 noise의 푸리에계수는 모두 제거한 것이다. 

남은 푸리에 계수만 가지고 푸리에 역변환을 하여 측정하고 자하는 원래의 파형을 얻을 수 있다.  

---