다운로드
작성자: admin 작성일시: 2017-08-01 16:10:43 조회수: 5025 다운로드: 135
카테고리: Python 태그목록:

푸리에 변환과 스펙트럼

음성(speech), 음악(music) 등의 음향(sound) 데이터에서 특징(feature)을 추출하는 방법인 푸리에 변환(Fourier transfrom)과 스펙트럼(spectrum)에 대해서 알아본다.

복소 지수함수

오일러 공식에 의해 지수부가 허수(imaginary number)인 복소 지수함수(complex exponential function)는 코사인 함수인 실수부와 사인 함수인 허수부의 합으로 나타난다.

$$ \exp(i\cdot x) = \cos x + i\sin x $$

다음은 주기가 $T$, 주파수가 $\frac{2\pi}{T}$인 복수 지수함수로서 주기가 $T$인 사인과 코사인의 조합이 된다.

$$ \exp \left( i\cdot 2\pi\frac{1}{T} t \right) = \cos \left(2\pi \dfrac{1}{T} t \right) +i\sin \left(2\pi \dfrac{1}{T} t \right) $$

따라서 주기가 $\frac{T}{n}$, 주파수가 $n\frac{2\pi}{T}$인 복수 지수함수는 다음과 같다.

$$ \exp \left( i\cdot 2\pi\frac{n}{T} t \right) = \cos \left(2\pi \dfrac{n}{T} t \right) +i\sin \left(2\pi \dfrac{n}{T} t \right) $$

푸리에 변환

주기 $T$를 가지고 반복되는(cyclic) 모든 함수 $f(t)$는 주파수와 진폭이 다른 몇 개의 사인 함수(정확히는 복수 지수함수)의 합으로 나타낼 수 있다. 이 사인 함수의 진폭을 구하는 과정을 푸리에 변환이라고 한다.

$$ y(t)=\sum_{k=-\infty}^\infty A_k \, \exp \left( i\cdot 2\pi\frac{k}{T} t \right) $$

이 식에서 $k$번째 사인 함수의 진폭 $A_n$은 다음 식으로 계산한다. 이것이 푸리에 변환이다.

$$ A_k = \frac{1}{T} \int_{-\frac{T}{2}}^\frac{T}{2} f(t) \, \exp \left( -i\cdot 2\pi \frac{k}{T} t \right) \, dt $$

이산 푸리에 변환

이산 푸리에 변환(Discrete Fourier Transform) 또는 DFT는 길이가 $N$인 이산시간 시계열 데이터

$$ y_0, y_1, \ldots, y_{N-1} $$

이 있을 때 이 이산시간 시계열이 주기 $N$으로 계속 반복된다고 가정하여 푸리에 변환을 한 것이다.

이 때 원래의 이산시간 시계열 데이터는 다음 주파수와 진폭이 다른 $N$개의 사인 함수의 합으로 나타난다.

$$ y_n = \frac{1}{N} \sum_{k=0}^{N-1} Y_k \cdot \exp \left( i\cdot 2\pi\frac{k}{N} n \right) $$

이 때 진폭 $Y_k$를 원래의 시계열 데이터에 대한 푸리에 변환값이라고 한다.

$$ Y_k = \sum_{n=0}^{N-1} y_n\cdot \exp \left( -i\cdot 2\pi\frac{k}{N} n \right) $$

FFT

FFT(Fast Fourier Transform)는 아주 적은 계산량으로 이산 푸리에 변환값을 계산하는 알고리즘을 말한다. 길이가 $2^N$인 시계열에만 적용할 수 있다는 단점이 있지만 보통의 DFT가 $O(N^2)$ 수준의 계산량을 요구하는데 반해 FFT는 $O(N\log_2 N)$ 계산량으로 DFT를 구할 수 있다.

예를 들어 다음과 같은 길이 $N=32$인 시계열 데이터를 푸리에 변환한다는 것은

In [1]:
N = 512
T = 1.0 / 44100.0
f1 = 697
f2 = 1209
t = np.linspace(0.0, N*T, N)
y1 = 1.1 * np.sin(2 * np.pi * f1 * t)
y2 = 0.9 * np.sin(2 * np.pi * f2 * t)
y = y1 + y2

plt.subplot(311)
plt.plot(t, y1)
plt.title(r"$1.1\cdot\sin(2\pi\cdot 697t)$")
plt.subplot(312)
plt.plot(t, y2)
plt.title(r"$0.9\cdot\sin(2\pi\cdot 1209t)$")
plt.subplot(313)
plt.plot(t, y)
plt.title(r"$1.1\cdot\sin(2\pi\cdot 697t) + 0.9\cdot\sin(2\pi\cdot 1209t)$")
plt.tight_layout()
plt.show()

실제로는 다음과 같이 계속 반복되는 시계열에 대해 푸리에 변환을 하는 것이다. 따라서 시계열의 시작 부분과 끝 부분이 너무 다르면 원래 시계열에는 없는 신호가 나올 수도 있는데 이를 깁스 현상(Gibbs phenomenon)이라고 한다.

In [2]:
y2 = np.hstack([y, y, y])
plt.plot(y2)
plt.axvspan(N, N * 2, alpha=0.3, color='green')
plt.xlim(400, 3 * N - 400)
plt.show()

scipy 패키지의 fftpack 서브패키지에서 제공하는 fft 명령으로 이 신호에 담겨진 주파수를 분석하면 다음과 같이 692HZ와 1211Hz 성분이 강하게 나타나는 것을 볼 수 있다. 이와 같은 플롯을 피리오도그램(periodogram)이라고 한다.

In [3]:
from scipy.fftpack import fft

yf = fft(y, N)
xf = np.linspace(0.0, 1.0/(2.0*T), N//2)

plt.subplot(211)
plt.plot(y)

plt.subplot(212)
plt.stem(xf, 2.0/N * np.abs(yf[0:N//2]))
plt.xticks(xf, rotation=40)
plt.xlim(0, 1989)

plt.show()

스펙트럼

보통 스펙트럼(spectrum)이라고 부르는 시계열 분석의 정확한 명칭은 파워 스펙트럼(power spectrum) 또는 스펙트럼 밀도(spectral density)이다.

푸리에 변환은 결정론적인 시계열 데이터를 주파수 영역으로 변환하는 것을 말하지만 스펙트럼(spectrum)은 확률론적인 확률과정(random process) 모형을 주파수 영역으로 변환하는 것을 말한다. 따라서 푸리에 변환과 달리 시계열의 위상(phase) 정보는 스펙트럼에 나타나지 않는다.

스펙트럼을 추정할 때 사용하는 방법 중의 하나전체 시계열을 짧은 구간으로 나눈 뒤 깁스 현상을 줄위기 위해 각 구간에 윈도우를 씌우고 FFT 계산으로 나온 값을 평균하는 방법이다. 보통은 로그 스케일로 표현한다.

In [4]:
import scipy.signal

f, P1 = sp.signal.periodogram(y, 44100, nfft=2**12)
f, P2 = sp.signal.welch(y, 44100, nfft=2**12)

plt.subplot(211)
plt.semilogy(f, P1)
plt.xlim(100, 1900)
plt.ylim(1e-8, 1e-1)

plt.subplot(212)
plt.semilogy(f, P2)
plt.xlim(100, 1900)
plt.ylim(1e-8, 1e-1)

plt.tight_layout()
plt.show()