푸리에 변환

시간이나 공간에 대한 함수를 시간 또는 공간 주파수 성분으로 분해하는 변환을 말한다. 종종 이 변환으로 나타난 주파수 영역에서 함수를 표현한 결과물을 가리키는 용어로도 사용된다.

: 임의의 입력 신호를 다양한 주파수를 갖는 주기함수들(sin, cos)의 합으로 분해하여 표현

ref: https://darkpgmr.tistory.com/171

$$f(x)=\int\limits_{-\infty}^{\infty}F(u)e^{j2\pi ux}du \rightarrow \text{푸리에 역변환 inverse Fourier transform}$$ $$F(u)=\int\limits_{-\infty}^{\infty}f(x)e^{-j2\pi ux}dx \rightarrow \text{푸리에 변환 Fourier transform}$$ $$j=\sqrt{-1}, f(x)=원본입력신호, e^{j2\pi ux}=\text{주파수 u인 주기함수성분}$$ $$F(u)= \text{해당 주기함수 성분의 계수(coefficient)}$$

$$\text{오일러 공식} e^{j\theta}=\cos\theta + j \sin\theta$$

$$e^{j2 \pi ux}=\cos 2\pi ux+j\sin 2\pi ux$$ $$\cos 2\pi ux, j\sin 2\pi ux \text{모두 주기period가 1/u, 주파수frequency가 u인 주기함수}$$ $$e^{j2 \pi ux}\text{= 주파수가 u인 정현파sinsusoidal wave의 복소지수함수 표현}$$ $$\text{정현파sinusoidal wave는 파형이 sin 또는 cos 함수인 파동wave}$$

  • 푸리에 역변환은 함수f(x)를 모든 가능한 주파수u의 주기함수들$e^{j2\pi ux}$의 일차결합으로 표현
  • 그 일차결합 계수 F(u)는 푸리에 변환으로 항상 주어질 수 있음
  • 푸리에 변환식을 푸리에 역변환에 대입했을떄 f(x)가 나옴을 증명해야 함.
    • 이 말은 임의의 신호 함수가 항상 주기함수들의 일차결합으로 분해될 수 있음을 증명하는 것
  • $e^{j2\pi ux}, u=0,\pm1, \pm2$는 모든 신호를 생성할 수 있는 직교 기저 함수orthogonal basis function
    • 여기서 u는 실수 전체 범위
    • 기저basis: 어떤 벡터 공간을 생성할 수 있는 일차 독립인 벡터들의 집합
    • 만일 기저 벡터들이 서로 수직orthogonal인 단위벡터라면 일차 결합 계수$a_i$는 dot product를 이용하여 $a_i=v\times v_i$로 계산
    • 어떤 벡터와 기저 벡터 dot product 시 이 벡터에 포함된 기저 성분의 계수가 얻어짐
  • 입력 신호 f(x)를 직교 기저 함수로 분해했을 떄의 F(u) 계수는 f(x)와 기저함수의 dot product 내적으로 계산될 수 있다.
  • 푸리에 변환은 f(x) 와 $e^{j2\pi ux}$의 함수 dot product이기 때문에 결과는 f(x)를 $e^{j2\pi ux}$들로 분해했을 때의 계수
    • 푸리에 변환에서 -가 붙은 이유는 복소수에서 dot product는 한 쪽에 conjugate 복소수를 취한 후 계산되기 때문

- 영상 신호에서 푸리에 변환

이미지는 2차원의 신호이기 때문에 2차원에서 정의되는 푸리에 변환 필요

$$f(x,y)=\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}F(u,v)e^{j2\pi (ux+vy)}dudv$$ $$F(u,v)=\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}f(x,y)e^{-2j\pi(ux+vy)}dxdy$$ $$F(u,v)\text{는 x축 방향으로 주파수 u, y축 방향으로 v인 주기함수 성분의 계수}$$

  • 이미지는 연속이 아닌 이산신호.
  • 한정된 유한 구간에서 정의되는 신호
  • 이산 데이터에 정의되는 푸리에 변환 필요

$$f(x,y)=\sum_{u=0}^{W-1}\sum_{v=0}^{H-1}F(u,v)e^{j2\pi (ux/W+vy/H)}\rightarrow\text{이산 데이터에서 푸리에 역변환}$$ $$F(u,v)=\frac{1}{WH}\sum_{x=0}^{W-1}\sum_{y=0}^{H-1}f(x,y)e^{-j2\pi(ux/W+vy/H)}\rightarrow\text{이산 데이터에서 푸리에 변환}$$ $$단, x=0,1,...,W-1 \dots u=0,1,...,W-1,\dots v=0,1,...H-1$$

  • 푸리에 역변환에서 $e^{j2\pi (ux/W+vy/H)}$는 x축 방향으로 주파수가 u/W, y축 방향으로 주파수가 v/H인 sinusoidal 주기 함수
  • 일반적인 푸리에 변환식과 달리 W와 H로 나눔 $\rightarrow$ 데이터가 정의된 구간을 하나의 단위 주리 unit period로 만드는 효과
    • 일종의 nomalization factor
  • 이미지를 신호로 해석하는 문제는 x,y축을 시간축으로 놓고 좌표의 변화에 따라 변하는 이미지 픽셀의 밝기 변화를 신호로 생각
  • 2D에서 정의되는 정현파의 모습은 모든 방향으로의 단면이 정현파가 되는 물결 형태의 파동 생각

이산 푸리에 변환에서 F(u,v)는 주파수 u,v의 성분이 아니라 주파수 u/W, v/H 성분에 대한 계수

W$\times$ H 인 이미지에서 이산 푸리에 변환의 F(u,v)는

  • x축 주파수가 u/W, y축 주파수가 w/H인 주기함수 성분에 대응
  • 주기로는 x축 방향 W/u 픽셀, y축 방향 H/v 픽셀인 주기성분을 나타냄
    • 주파수 공간에서 특정 F(u,v)값이 높게 나타났다면 원래의 이미지 공가에서는 x축 방향으로의 주기가 W/u 픽셀, y축 방향 주기가 H/v 픽셀인 주기성 성분이 존재

참고 $$f(x)=\sum_{u=0}^{W-1}F(u)e^{j2\pi ux/W},ㅌ=0,1,\dots ,W-1\rightarrow\text{1차원 이산 데이터에서 푸리에 역변환}$$ $$F(u)=\frac{1}{WH}\sum_{x=0}^{W-1}f(x)e^{-j2\pi ux/W},u=0,1,\dots,W-1\rightarrow\text{이산 데이터에서 푸리에 변환}$$ -1차원에서 이산 푸리에 변환 식은f(x)가 주기함수 일때만 성립

- 푸리에 스펙트럼과 페이즈

$$F(u,v)=R(u,v_+iI(u,v)$$

  • 푸리에 변환으로 구한 F(u,v)는 복소수complex number이고, 실수부real과 허수부Imaginary로 구성
  • 복소수 F(u,v)의 크기 |F(u,v)|는 푸리에 변환의 스펙트럼spectrum 또는 magnitude 라 부름
  • F(u,v)의 각도$\phi$ 를 페이즈phase angle 또는 phase spectrum이라 부름 $$|F(u,v)|=[R^2(u,v)+I^2(u,v)]^{1/2}$$ $$\phi(u,v)=\tan^{-1}\big[\frac{I(u,v)}{R(u,v)}\big]$$

푸리에 스펙트럼 Fourier Spectrum

  • 해당 주파수 성분이 원 신호에 얼마나 강하게 포함되어 있는지
  • 푸리에 변환 및 역변환으로 얻은 W$\times$H의 F(u,v), |F(u,v)|를 픽셀값으로 잡으면 푸리에 스펙럼을 원본 이미지와 동일한 크기의 이미지로 시각화

문제점

  • 푸리에 스펙트럼의 저주파 영역은 매우 큰 값을 갖지만 다른 영역은 0에 가까운 값을 갖기 떄문에 푸리에 스펙트럼을 그대로 이미지로 시각화하면 검은 바탕 위에 흰점 하나만 존재하는 형태
    • 해결: 스펙트럼을 이미지로 표현할 때는 스펙트럼에 log를 취하는 것이 일반적
  • 원래 스펙트럼 이미지는 모서리로 갈수록 값이 높아지기 때문에 스펙트럼 형태를 파악하기 힘듦
    • 해결: 원점이 중심에 오도록 스펙트럼 위치를 이동shift시킨 형태의 이미지를 사용하는 것이 일반적$\rightarrow$ 푸리에 스펙트럼 이미지라 함은 shift된 스펙트럼 이미지를 가리킴
      • 푸리에 스펙트럼이 원점 대칭인 주기함수이기 때문에 가능한 shift!

푸리에 변환의 페이즈phase

  • 페이즈phase: 파형의 시점이 어디인지?
    • sin파와 cos파는 90도의 phase차이가 존재하는 동일한 파형으로 볼 수 있음
  • 푸리에 변환의 관점에서 보면 phase는 원본 신호를 주기 신호로 분해했을때 각 주기 성분의 시점이 어디인지(즉, 각 주기 성분들이 어떻게 줄을 맞춰서 원본 신호를 생성했는지)를 나타내는 요소가 된다.

$$F(u) = R(u)jI(u)=|F(u)|e^{j\phi (u)}$$ $$1차원의 경우임$$

$$f(x) = \int\limits_{-\infty}^{\infty}|F(u)|e^{j(2\pi u x + \phi (u))}du$$

푸리에 변환 식에 대입했을때 phase term이 주기함수 성분의 시점을 조절하는 term이 된다.

  • 즉, 푸리에 계수 F(u)에 대응되는 주기함수 성분의 강도를 나타내는 스펙트럼 정보 |F(u)|와 시점을 조절하는 phase의 정보$\phi (u)$ 가 함께 포함되어 있음을 알 수 있다.

ref: https://datascienceschool.net/03%20machine%20learning/03.03.02%20%ED%91%B8%EB%A6%AC%EC%97%90%20%EB%B3%80%ED%99%98%EA%B3%BC%20%EC%8A%A4%ED%8E%99%ED%8A%B8%EB%9F%BC.html?highlight=%ED%91%B8%EB%A6%AC%EC%97%90

정현파 조합

모든 신호는 주파수frequency와 크기magnitude, 위상phase이 다른 정현파simusolida signal의 조합으로 나타낼 수 있다. 푸리에 변환은 조합된 정현파의 합 신호에서 그 신호를 구성하는 정현파들을 각각 분리해내는 방법이다.

import numpy as np
import matplotlib.pyplot as plt


N = 1024
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()

복소 지수함수

오일러 공식에 의해 지수부가 허수imaginary number인 복소 지수함수complex exponential function는 코사인 함수인 실수부와 사인 함수인 허수부의 합으로 나타난다. $$exp(i \cdot x) = \cos z + i \sin x$$

다음은 주기가 T, 주파수가 $\frac{2\pi}{T}$인 복소 지수함수로서 주기가 T인 sin과 cos의 조합이 된다. $$exp\big(i \cdot 2\pi \frac{1}{T}t \big)=\cos \big(2\pi frac{1}{T}t \big) + i\sin\big(2\pi frac{1}{T}t \big)$$

따라서 주기가 $frac{T}{n}$, 주파수가 $n\frac{2\pi}{T}$인 복소 지수함수는 다음과 같다. $$exp\big(i \cdot 2\pi \frac{n}{T}t \big)=\cos \big(2\pi frac{n}{T}t \big) + i\sin\big(2\pi frac{n}{T}t \big)$$

푸리에 변환

주기 T를 가지고 반복되는cycle 모든 함수 $y(t)$는 주파수와 진폭이 다른 몇 개의 sin 함수(정확히는 복소 지수함수)의 합으로 타나낼 수 있다. 이 sin 함수의 직폭을 구하는 과정을 푸리에 변환Fourier Transform이라고 한다. $$y(t)=\sum^{\infty}_{k=-\infty} A_k exp\big( i \cdot 2\pi frac{k}{T}t \big)$$

이 식에서 k번째 함수인 sin함수의 진폭 $A_n$은 다음 식으로 계산한다. 이것이 푸리에 변환이다. $$A_n = \frac{1}{T}\int\limits_{-\frac{T}{2}}^{\frac{T}{2}}y(t)exp\big( -i \cdot 2\pi \frac{k}{T}t \big)$$

이산 푸리에 변환

이산 푸리에 변환Discrete Fourier Transform 또는 DFT는 길이가 N인 이산시간 시계열 데이터 $y_0,y_1, \dots , y_{N-1}$이 있을때, 이 이산시간 시계열이 주기 N으로 계속 반복된다고 가정하여 푸리에 변환을 한 것이다.

이 때 원래의 이산시간 시계열 데이터는 다음 주파수와 진폭이 다른 N개의 sin 함수의 합으로 나타낸다. $$y_n=\frac{1}{N} \sum^{N-1}_{k=0}Y_k \cdot exp \big( -i \cdot 2\pi \frac{k}{N}n \big)$$

이 때 진폭 Y_k를 원래의 시계열 데이터에 대한 푸리에 변환값이라고 한다. $$Y_k = \sum^{N-1}_{n=0}y_n \cdot exp \big( i \cdot 2\pi \frac{k}{N}n \big)$$

고속 푸리에 변환

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

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

y2 = np.hstack([y, y, y])

plt.subplot(211)
plt.plot(y2)
plt.axvspan(N, N * 2, alpha=0.3, color='green')
plt.xlim(0, 3 * N)

plt.subplot(212)
plt.plot(y2)
plt.axvspan(N, N * 2, alpha=0.3, color='green')
plt.xlim(900, 1270)

plt.show()

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

from scipy.fftpack import fft

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

plt.stem(xf, 2.0/N * np.abs(yf[0:N//2]))
plt.xlim(0, 3000)

plt.show()

DCT

DCT Dicrete Cosine Tranform은 DFT와 유사하지만 기저함수로 복소 지수함수가 아닌 cos 함수를 이용한다. DFT보다 계산이 간단하고 실수만 출력한다는 장점이 있어서 DFT 대용으로 많이 사용한다. $$Y_k=\sum^{N-1}_{n=0} y_n \cdot \cos \big( 2\pi frac{k}{N}\big( \frac{2n+1}{4} \big) \big)$$ 

from scipy.fftpack import dct

dct_type = 2
yf2 = dct(y, dct_type, N)

plt.subplot(311)
plt.stem(np.real(yf))
plt.title("DFT real") # 실수부

plt.subplot(312)
plt.stem(np.imag(yf))
plt.title("DFT imaginary") # 허수부

plt.subplot(313)
plt.stem(np.abs(yf2))
plt.title("DCT")

plt.tight_layout()
plt.show()

스펙트럼

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

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

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

import scipy.fft as fft

import scipy.signal

f, P = scipy.signal.periodogram(y, 44100, nfft=2**12)

plt.subplot(211)
plt.plot(f, P)
plt.xlim(100, 1900)
plt.title("Linear Scale") # 선형 스케일

plt.subplot(212)
plt.semilogy(f, P)
plt.xlim(100, 1900)
plt.ylim(1e-5, 1e-1)
plt.title("Log Scale") # 로그 스케일

plt.tight_layout()
plt.show()

STFT

STFT Short-Time Fourier Transform 는 주파수 특성이 시간에 따라 달라지는 사운드를 분석하기 위한 방법이다. 시계열을 일정한 시간 구간으로 나누고 각 구간에 대해 스펙트럼을 구한 데이터이다. 시간-주파수의 2차원 데이터로 나타낸다.

librosa package

파이썬으로 STFT 스펙트럽 분석을 하려면 librosa package를 사용한다.

주피터 노트북에서 librosa 패키지를 사용할 때는 jupyter_notebook_config.py 파일의 iopub_data_rate_limit 설정을 10000000 정도로 크게 해야 한다.

import librosa
import librosa.display

D = np.abs(librosa.stft(y))
librosa.display.specshow(librosa.amplitude_to_db(D, ref=np.max), y_axis='linear', x_axis='time')
plt.title('Dual Tone')
plt.ylim(0, 4000)
plt.show()

/home/csy/anaconda3/envs/csy/lib/python3.8/site-packages/librosa/core/spectrum.py:222: UserWarning: n_fft=2048 is too small for input signal of length=1024
  warnings.warn(

sr_octave, y_octave = scipy.io.wavfile.read("octave.wav")

D_octave = np.abs(librosa.stft(y_octave))
librosa.display.specshow(librosa.amplitude_to_db(D_octave, ref=np.max), sr=sr_octave, y_axis='linear', x_axis='time')
plt.title('Octave')
plt.ylim(0, 2000)

ㅡ음악 파일 없음..

멜 스펙트럼

멜 스펙트럼은 주파수의 단위를 다음 공식에 따라 멜 단위mel unit으로 바꾼 스펙트럼을 말한다. $$m=2595 log_10 \big( 1+ \frac{f}{700} \big)$$ 

S_octave = librosa.feature.melspectrogram(y=y_octave, sr=sr_octave, n_mels=128)
librosa.display.specshow(librosa.power_to_db(S_octave, ref=np.max), sr=sr_octave, y_axis='mel', x_axis='time')
plt.ylim(0, 4000)
plt.show()

MFCC

![]("https://images1.programmersought.com/188/d5/d590aa0094ee101ec8ad4e811f06e1fc.JPEG")

/bin/bash: -c: line 0: syntax error near unexpected token `"https://images1.programmersought.com/188/d5/d590aa0094ee101ec8ad4e811f06e1fc.JPEG"'
/bin/bash: -c: line 0: `[]("https://images1.programmersought.com/188/d5/d590aa0094ee101ec8ad4e811f06e1fc.JPEG")'

MFCC Mel-frequency cepstral coefficients는 Mel Scale Spectrum을 40개의 주파수 구역band으로 묶은 뒤에 이를 다시 푸리에 변환하여 얻은 계수이다. 스펙트럼이 어떤 모양으로 되어 있는지를 나타내는 특성값이라고 생각할 수 있다.

y, sr = librosa.load(librosa.util.example_audio_file(), offset=30, duration=5)

plt.plot(y[1000:5000])
plt.show()

from IPython.display import Audio, display
Audio(y, rate=sr)

S = librosa.feature.melspectrogram(y=y, sr=sr, n_mels=128, fmax=8000)
librosa.display.specshow(librosa.power_to_db(S, ref=np.max), sr=sr, y_axis='mel', x_axis='time')
plt.show()

mfccs = librosa.feature.mfcc(y=y, sr=sr, n_mfcc=40)
librosa.display.specshow(mfccs, x_axis='time')
plt.title('MFCC')
plt.tight_layout()
plt.show()

푸리에 변환; 단순한 파형들의 합성이 어떤 주파수나 세기로 성립되어 있는지를 수학적으로 구하는 방법

  • 파형을 이용해 스펙트럼을 구하는 방법
  • 푸리에 역변환; 스펙트럼으로 파형구하기..

직교; 원주 위에서 위치관계가 90도(직각)만큼 차이가 나는 상태

  • $\sin mx$와 $\cos mx$ 는 직교하고 있다

푸리에 급수

$$F(x) = \frac{1}{2}a_0 + \sum^{\infty}_{n=1} (a_n \cos nx + b_n \sin nx)$$

푸리에 계수

$$a_n = \frac{1}{\pi} \int^{2\pi}_{0} F(x) \cos nx dx$$

$$b_n = \frac{1}{\pi} \int^{2\pi}_{0} F(x) \sin nx dx$$

$$a_0 = \frac{1}{2\pi} \int^{2\pi}_{0} F(x) dx$$

  • (cos 0x = 1) 곱해져서 안 보여~