Chapter4 Vectorization: Notebook6 2-dim Vectorization Example3

NumPy Master Class¶

Chapter4 Vectorization¶

Notebook6 2-dim Vectorization Example3¶

이번 2차원 Vectorization에서는 notebook3에서 다뤘던 DTFT를 하나의 주파수가 아닌 주파수 영역에서 구해보려고 한다. 먼저 Foureir Transform의 식을 임의의 주파수에서 구하는 식을 보면 다음과 같다.

$X(f) = \sum_{n= -\infty}^{\infty} x[n]*e^{-j2 \pi fn}$

그러면 우리의 신호를 다음과 같이 만들어보자.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import time

from scipy.io import wavfile
fs, data = wavfile.read('./test_audio.wav')
data = (data[:,0] + data[:,1])/2
fig, ax = plt.subplots(figsize = (15,6))
ax.plot(data)

/home/shinks/anaconda3/envs/pytorch/lib/python3.6/site-packages/ipykernel_launcher.py:6: WavFileWarning: Chunk (non-data) not understood, skipping it.
  

Out[1]:

[<matplotlib.lines.Line2D at 0x7f3fd85919e8>]

이 audio에 대하여 -20Hz부터 20Hz까지의 DTFT를 2중 for loop을 이용하여 구하고, frequency response는 하나의 for loop을 이용해 구해보자.

In [2]:

n_f, n_s = 100, data.shape[0]

f_range = np.linspace(-20, 20, n_f).reshape(n_f,1)
Xf = np.zeros(shape = (f_range.shape), dtype = np.complex)

In [4]:

tic = time.time()
PI = np.pi

for f_idx in range(n_f):
    f = f_range[f_idx]
    for s_idx in range(n_s):
        s = data[s_idx]
        Xf[f_idx] += s * np.exp(-2j*PI*f*s_idx)
        
for f_idx in range(n_f):
    Xf[f_idx] = np.abs(Xf[f_idx])
    
toc = time.time()
for_time = toc - tic
print("Elapsed Time:", for_time, 'sec')

plt.figure(figsize = (10,5))
plt.plot(Xf)
plt.yscale('log')

Elapsed Time: 152.4426200389862 sec

/home/shinks/anaconda3/envs/pytorch/lib/python3.6/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

그러면 위의 과정을 Vectorization으로 구해보자.

다음 연산은 여러 단계로 나누는 것이 더 알아보기 좋으므로 Step으로 나눠서 알아보자.

$X(f) = \sum_{n= -\infty}^{\infty} x[n]*e^{-j2 \pi fn}$

전체 단계 중에서 다음과 같이 나눌 수 있다.

먼저 가장 크게 Vectorized operation이 진행되는 부분은 전체 f와 data를 곱하는 부분인 step1이다.

step1: $fn$

step2: $-j2 \pi fn$

step3: $e^{-j2 \pi fn}$

step4: $x[n]*e^{-j2 \pi fn}$

step5: $\sum_{n= -\infty}^{\infty} x[n]*e^{-j2 \pi fn}$

In [5]:

tiled_data = np.tile(data, (n_f, 1))
print("tiled_data.shape:", tiled_data.shape) # (f 개수, data 개수)
print("f_range.shape:", f_range.shape)

step1 = tiled_data*f_range
print("step1.shape:", step1.shape)

tiled_data.shape: (100, 268237)
f_range.shape: (100, 1)
step1.shape: (100, 268237)

위와 같이 broadcasting을 이용하여 f의 값들을 차례대로 data들에 곱해준다.

In [6]:

step2 = -2j*PI*step1
print("step2.shape:", step2.shape)

step2.shape: (100, 268237)

In [7]:

step3 = np.exp(step2)
print("step3.shape:", step3.shape)

step3.shape: (100, 268237)

In [8]:

print("step3.shape:", step3.shape)
print("data.shape:", data.shape)

step4 = step3*data
print("step4.shape:", step4.shape)

step3.shape: (100, 268237)
data.shape: (268237,)
step4.shape: (100, 268237)

In [9]:

step5 = np.abs(np.sum(step4, axis = 1))
print(step5.shape)

(100,)

위의 과정에서 frequency response를 구하기 위해 절대값까지 같이 취해준다.

위의 과정을 모두 종합하면 다음과 같다.

In [10]:

tic = time.time()
PI = np.pi

DTFT = np.abs(np.sum((np.exp(-2j*PI*(np.tile(data, (n_f, 1)))*f_range))*data, axis = 1))

toc = time.time()
vec_time = toc - tic
print("Elapsed Time:", vec_time, 'sec')

plt.figure(figsize = (10,5))
plt.plot(DTFT)
plt.yscale('log')

Elapsed Time: 1.9441606998443604 sec

In [11]:

print("Elapsed Time Ratio:", for_time / vec_time)

Elapsed Time Ratio: 78.41050384939372

즉, 하나의 test audio file에 대해 FFT를 할 때도 78배의 연산속도 차이가 났다.

위의 data들의 shape을 살펴보면 다음과 같다.

In [12]:

print("data.shape:", data.shape)
print("f_range.shape:", f_range.shape)

data.shape: (268237,)
f_range.shape: (100, 1)

위와 같이 하나의 30초짜리 audio file도 26만 point가 나오므로 실제로 우리가 다루게 되는 data는 훨씬 더 클 것을 예상할 수 있다.

따라서 vectorization이 더더욱 중요해질 것은 분명해 보인다.