[기계진동] Spectral Leakage와 Aliasing

이번 포스팅에서는 Spectral Leakage와 Aliasing이 무엇인지에 대해

 

Python으로 직접 해보면서 알아보도록 합시다!

 

 Spectral Leakage

Spectral Leakage란, 신호를 주파수 분석 했을 때

 

원래의 신호에는 포함되어 있지 않은 주파수 성분이 관측되는 현상을 말합니다.

 

원인은 신호의 주파수 성분이 샘플링 주파수의 n배와 매칭되지 않을 때 나타나죠!

 

예제를 통해 알아볼까요? 아래 함수를 FFT한다고 해봅시다

 

함수 f(x)는 55Hz, 120Hz 성분의 사인함수를 포함하고 있기 때문에

 

FFT했을 때 기대되는 모양은 55Hz, 120Hz의 주파수에서 peak가 형성될 것입니다.

 

하지만 만약 t=1/400,  N=40,  T=1/10의 조건으로 FFT하게되면 어떨까요?

(t: sampling time, N: number of data, T=period)

 

그림1. FFT 조건: t=1/400,    N=40,    T=1/10

 

주파수 성분의 peak 주변으로 원 함수에 포함되어있지 않은 주파수 성분들이 나타납니다.

 

특히 55Hz 주변이 두드러지네요.

 

1/T=10Hz이기 때문에, 55Hz는 10의 배수가 아니기 때문입니다.

 

그럼 이를 방지하려면? 주기의 역수를 55Hz의 배수로 만들어주면 되니까

 

데이터의 갯수를 늘리면 되겠네요. 데이터 갯수를 늘려볼까요?

 

N=80으로 설정해보겠습니다.

 

그림2. FFT 조건: t=1/400,    N=80,    T=1/5

그림 2와 같이 깔끔하게 나타나는 그래프를 확인할 수 있습니다.

 

 

 Aliasing

Aliasing은 입력신호의 최대주파수보다 Nyquist frequency가 적어서,

 

신호를 복원하는 과정에서 최대 주파수 성분을 구분하지 못하는 현상을 말합니다.

 

위의 예제 함수를 통해 알아보죠.

 

t=1/200,  N=40,  T=1/5의 조건으로 FFT를 한다고 해봅시다.

 

이 때, Nyquist frequency fn = 1/(2t) = 100Hz 입니다.

그림 3. FFT 조건: t=1/200,    N=40,    T=1/5

원래의 신호는 55Hz 성분과 120Hz 성분을 포함하고 있었는데

 

FFT 결과는 어찌된 일인지 120Hz는 온데간데없고 80Hz에 peak가 나타났습니다.

 

원래의 120Hz 신호가 Nyquist frequency보다 높았기 때문에 잡아내지 못한거죠!

 

만약 이 신호를 다시 시간역으로 복원한다면?

네... 원래의 하늘색 신호와 매우 다른 모습을 볼 수 있습니다.

 

실전에서 만약 이렇게 되면 큰일나겠죠?

 

무작정 sampling frequency를 키우기만 하는 것은 경제적이지 못하지만,

 

적절한 필터 사용을 통해 경제적인 Nyquist frequency를 선정해야 하겠습니다.

 

 

 

* 본문의 그래프는 MATLAB을 통해 따온 그래프입니다.

  아래 파이썬 코드를 작동했을 때도 동일한 결과가 나와야 할 것 같은데...

  유효숫자때문인 것 같기도 하고 MATLAB과 같은 결과가 나오질 않네요 ㅠㅠ

  아시는 분은 댓글로 도움부탁드립니다...

 

import numpy as np
import matplotlib.pyplot as plt

dt = 1/200
Fs = 1/dt
N = 40
t_end = N*dt
len_data = int(t_end/dt + 1)
t = np.arange(len_data)*dt
n = np.arange(len_data)
T = len_data/Fs
freq = n/T

sin_120Hz = np.array([])
sin_55Hz = np.array([])
fn_sin = np.array([])

for i in range(0, len_data):
    sin_120Hz = np.append(sin_120Hz, np.sin(2.0 * np.pi * 120.0 * dt * i))
    sin_55Hz = np.append(sin_55Hz, np.sin(2.0 * np.pi * 55.0 * dt * i))

fn_sin = sin_55Hz + sin_120Hz

Y_raw = np.fft.fft(fn_sin)/len_data

y = np.fft.ifft(Y_raw)*len_data

plt.figure()
plt.plot(t, fn_sin,'b')
plt.axis([0, t_end, -2, 2])
plt.grid(True)
plt.xlabel('Second [s]')
plt.ylabel('Amplitude [m]')
plt.plot(t, y,'r--')
plt.axis([0, t_end, -2, 2])
plt.grid(True)
plt.xlabel('Second [s]')
plt.ylabel('Amplitude [m]')

plt.figure()
plt.stem(freq, np.abs(Y_raw)*2)
plt.axis([0, Fs, 0, 2])
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [m]')
plt.grid(True)

plt.show()