Fourier Transform

2. Fourier Transforms

푸리에 변환(Fourier transform)을 직관적으로 설명하면 푸리에 변환은 임의의 입력 신호를 다양한 주파수를 갖는 주기함수(복수 지수함수)들의 합으로 분해하여 표현하는 것이다. 그리고 각 주기함수들의 진폭을 구하는 과정을 퓨리에 변환이라고 한다.

 

2-1. Fourier Series

1807년 푸리에는 임의의 함수를 삼각함수의 (무한) 선형조합으로 표현할 수 있다고 주장하였다.(참고)

 

힐버트 공간 F=L2([0,2π])상에서 지수함수계 E={2π​1​ejnt:n=0,±1,±2,⋯​}는 완비정규직교 기저집합이므로 [0,2π]에서 정의되는 주기함수들은 다음과 같이 푸리에 급수로 표현된다.

x(t)=n=−∞∑∞​an​2π​1​ejntan​=⟨f,pn​⟩=∫02π​x(t)2π​1​e−jntdt

먼저 첫번째 식을 보면n는−∞~∞의 범위를 가지고 움직인다. 즉, 어떠한 신호가 서로 다른 주기함수들의 합으로 표현되는데, 그 주기함수는 무한개가 있어야 함을 의미한다.
한가지 더 주목할 점은 지수함수계의 함수들이 모두 직교한다는 점이다. 즉,x(t)는 기저를 E로 두면,(a−∞​,⋯,a∞​)라는 좌표로 표현할 수 있게 된다.

2-2. DFT (Discrete Fourier Transform)

□ Fourier Transform

위 결과를 주기가 무한대인 일반적인 함수로 일반화하면, Fourier Transform 식을 얻을 수 있다.

InverseFourierTransform : x(t)=2π1​∫−∞∞​X(jω)ejωtdω

FourierTransform : X(jω)=∫−∞∞​x(t)e−jωtdt

즉, time domain에서 연속적이고 무한한 길이의 신호를 frequency domain에서 연속적이고 무한한 주파수로 나타낼 수 있다.

 

□ Discrete Fourier Transform

우리가 sampling한 신호는 바로 시간의 간격과 소리의 amplitude가 모두 discrete한 데이터이다. 따라서 위의 푸리에 변환 식을 Discrete한 영역으로 바꾸어야 한다. 연속적인 신호에서 N개 샘플링한 이산적인 신호에서는 다음과 같이 푸리에 급수로 표현된다.

InverseDiscreteFourierTransform : x[t]=N1​k=0∑N−1​X[k]ejωk​n

DiscreteFourierTransform : X[k]=n=0∑N−1​x[n]e−jωk​n

우리가 가진 신호x[t]에서, 이산 시계열 데이터가 주기N으로 반복한다고 할때, DFT는 주파수와 진폭이 서로 다른 N개의 사인 함수의 합으로 표현이 가능하다.

위 식을 보면 k의 range가 0부터 N−1로 변화했음을 알 수 있다. 이때 Spectrum X[k]를 원래의 시계열 데이터에 대한 퓨리에 변환값이라 한다.

 

 

□ Zero Padding

  첫 번째 행의 입력 신호의 길이가 짧기 때문에 DFT로 나온 스펙트럼 또한 적은 개수의 주파수만을 표현할 수 있다. 즉, 신호의 길이가 짧아서 DFT 된 신호가 정확히 표현되지 않은 경우이다. 좀 더 정확한 주파수 영역의 신호를 얻기 위해서는, 위의 예시와 같은 문제점을 해결해야한다. 이 때 사용하는 방법이 바로 zero-padding 이다. zero padding은 x(n)의 뒤에 0을 여러개 덧붙여서 강제로 x(n)의 길이를 늘려줌으로써 DFT된 X(m)을 더 촘촘하게 나타내는 방법이다.

 

  신호의 뒤에 0을 덧붙여 DFT 를 수행한다면, 샘플 수가 늘어난 X(m)을 얻을 수 있다. 이러한 방식으로 고밀도의 스펙트럼을 얻을 수 있다 신호의 샘플 수는 늘어나게 되어 X(m)이 더욱 촘촘한 간격의 한주기 신호로 생성되었고, 단순히 0을 덧붙였으므로 신호의 특성이나 분해능에는 영향을 미치지 않는다. 따라서, zero padding은 고밀도의 스펙트럼을 얻게 해주지만 고분해능의 스펙트럼을 얻을 수 있는 것은 아니다.

 

 

위 그림과 같이 실제로 신호의 길이가 길어진 경우와 단순히 zero padding을 붙여서 신호의 길이를 길게 만든 경우의 스펙트럼은 다르게 나타난다. 신호의 길이가 길어진 경우 main lobe가 점점 좁아지는 것을 확인할 수 있다.

신호의 길이가 길어지면 원래 신호의 주파수가 어떻게 분포하는지를 정확하게 파악할 수 있다. 반면에 zero padding을 붙인 경우는 위에서 볼 수 있듯이 고밀도로 간격은 좁아졌지만 main lobe의 폭은 변화가 없기 때문에 고분해능이라 보기는 힘들다.

Code: DFT

import numpy as np
import matplotlib.pyplot as plt

def DFT(x):
  N = len(x)
  X = np.array([])
  nv = np.arange(N)
  for k in range(N):
    s = np.exp(1j*2*np.pi*k/N*nv)
    X = np.append(X, sum(x*np.conjugate(s)))
  return X

def plot_magnitude(DFT, ax, N):
  magnitude = np.abs(DFT)
  freq_axis = 2 * np.pi * np.arange(magnitude.shape[-1]) / N
  ax.stem(freq_axis, magnitude, 'k', markerfmt = 'ko', basefmt = " ")
  ax.set(title="Magnitude of DFT", xlim = (0, np.pi), xticks=[0, np.pi / 2, np.pi], xticklabels=[r'$0$', r'$\pi/2$', r'$\pi$'], xlabel = 'Frequency (rad)', ylabel='Manigude')

def plot_phase(DFT, ax, N):
  phase = np.angle(DFT)
  freq_axis = 2 * np.pi * np.arange(DFT.shape[-1]) / N
  ax.stem(freq_axis, phase, 'k', markerfmt = 'ko', basefmt = " ")
  ax.set(title="Phase of DFT", xlim = (0, np.pi), xticks=[0, np.pi / 2, np.pi], xticklabels=[r'$0$', r'$\pi/2$', r'$\pi$'], xlabel = 'Frequency (rad)', 
  ylabel='Phase (rad)', ylim=(-np.pi, np.pi), yticks=[-np.pi, 0, np.pi], yticklabels=[r'$-\pi$', r'$0$', r'$\pi$'])

num_samples = 47
sinusoid = np.cos(np.pi * 2 / 5.5 * np.arange(num_samples))
X = DFT(sinusoid)
fig, axes = plt.subplots(ncols=3)
axes[0].stem(np.arange(num_samples), sinusoid, 'k', markerfmt = 'ko', basefmt = " ")
axes[0].set(title="Discrete signal", xlabel="Discrete index", ylabel="Amplitude")
plot_magnitude(X, axes[1], sinusoid.shape[-1])
plot_phase(X, axes[2], sinusoid.shape[-1])
fig.set_size_inches(18, 6)
fig.tight_layout()
plt.show()

 

DFT의 실제 구현은 FFT를 이용하기 때문에 O(nlogn)의 시간복잡도를 가진다.(참고)

 

2-3. STFT (Short Time Fourier Transform)

spectrogram vs spectrum

import librosa
import IPython.display as ipd

wav, sr = librosa.load(librosa.ex("brahms"), duration = 3)
ipd.display(ipd.Audio(wav, rate=sr))
S = librosa.core.stft(wav, n_fft=1024, hop_length=512, win_length=1024)
log_S = librosa.power_to_db(np.abs(S)**2, ref=np.max)

fig, axes = plt.subplots(ncols=2)
plot_magnitude(np.fft.rfft(wav), axes[0], wav.shape[-1])

librosa.display.specshow(log_S, sr=sr, x_axis='time', y_axis='hz', ax = axes[1])
axes[1].set_title("STFT")
fig.set_size_inches(12, 4)
fig.tight_layout()
plt.show()

 

 왼쪽 그림이 Spectrum이고, 오른쪽 그림이 Spectrogram이다. 둘의 가장 큰 차이는 X축이 Time이냐, Frequency냐의 차이이다.

 Spectrum은 시간 축이 없기 때문에 특정 순간의 소리의 에너지를 분석한 것이다. 반면에 spectrogram은 시간에 따른 소리의 변화를 시각화한 것이다.

 대부분의 신호는 시간에 따라 주파수가 변하게 된다. 그러나 FFT는 어느 시간대에 주파수가 변하는지 알 수 없다. 즉 FFT를 하면 Time domain에 대한 정보가 사라진다. 이러한 한계를 극복하기 위해서, STFT는 시간을 프레임별로 나눠서 FFT를 수행한다.

 

□ STFT

X(l,k)=n=0∑N−1​w(n)x(n+lH)expN−2πkn​

STFT는 주파수의 특성이 시간에 따라 달라지는 사운드를 분석하는 방법이다. speech signal의 특성이 상대적으로 느리게 변한다고 가정하고, 짧은 시간간격 안에서는 stationary하다고 가정한다. 짧은 시간간격으로 나누는 framing과정을 거쳐서 short-time에 대한 processing을 하게된다. 시계열 데이터를 일정한 시간 구간 (window size)로 나누고, 각 구간에 대해서 스펙트럼을 구하는 데이터이다. 따라서 STFT를 수행하면, (time, frequncy, magnitude) 꼴로 spectrogram이 나온다. 시각화할 때는 일반적으로 시간을 x축으로, 주파수를 y축으로, 크기를 색깔로 표현한다.

 

□ Parameter of STFT

 ·  N : FFT size

   - Window를 얼마나 많은 주파수 밴드로 나누는가

 ·  w(n) : Window function

   - 일반적으로 Hann window가 쓰인다.

·  n : Window size

   - Window 함수에 들어가는 Sample의 양이다.

   - 작을수록 Low-frequency resolution을 가지게 되고, high-time resolution을 가진다.

   - 길수록 High-frequency, low time resolution을 가진다.

 ·  H : Hop size

   - 윈도우가 겹치는 사이즈이다. 일반적으로는 1/2 혹은 1/4정도를 겹치게 한다.