scipy.fftpack 모듈은 fast Fourier transforms (FFT)을 계산하고 이를 처리 할 수 있는 유틸리티를 제공합니다. 주요 기능은 다음과 같습니다.
scipy.fftpack.fft ()를 사용하여 FFT를 계산
scipy.fftpack.fftfreq ()를 사용하여 샘플링 주파수 생성
scipy.fftpack.ifft ()는 주파수 공간에서 신호 공간으로 inverse FFT를 계산
다음 예제를 실습해 보겠습니다. 이 예제는 signal의 FFT의 힘을 플롯하고 inverse FFT를 사용하여 signal을 재구성하는 것입니다. 이 예제는 scipy.fftpack.fft(), scipy.fftpack.fftfreq() 및 scipy.fftpack.ifft()를 보여줍니다.
import numpy as np
from scipy import fftpack
from matplotlib import pyplot as plt
#신호를 생성합니다.
# Seed the random number generator
np.random.seed(1234)
time_step = 0.02
period = 5.
time_vec = np.arange(0, 20, time_step)
sig = (np.sin(2 * np.pi / period * time_vec) + 0.5 * np.random.randn(time_vec.size))
plt.figure(figsize=(6, 5))
plt.plot(time_vec, sig, label='Original signal')
#FFT의 Power를 계산합니다.
# The FFT of the signal
sig_fft = fftpack.fft(sig)
# And the power (sig_fft is of complex dtype)
power = np.abs(sig_fft)
# The corresponding frequencies
sample_freq = fftpack.fftfreq(sig.size, d=time_step)
# Plot the FFT power
plt.figure(figsize=(6, 5))
plt.plot(sample_freq, power)
plt.xlabel('Frequency [Hz]')
plt.ylabel('plower')
# Find the peak frequency: we can focus on only the positive frequencies
pos_mask = np.where(sample_freq > 0)
freqs = sample_freq[pos_mask]
peak_freq = freqs[power[pos_mask].argmax()]
# Check that it does indeed correspond to the frequency that we generate
# the signal with
np.allclose(peak_freq, 1./period)
# An inner plot to show the peak frequency
axes = plt.axes([0.55, 0.3, 0.3, 0.5])
plt.title('Peak frequency')
plt.plot(freqs[:8], power[:8])
plt.setp(axes, yticks=[])
# scipy.signal.find_peaks_cwt can also be used for more advanced peak detection
#모든 high frequencies를 제거합니다.
high_freq_fft = sig_fft.copy()
high_freq_fft[np.abs(sample_freq) > peak_freq] = 0
filtered_sig = fftpack.ifft(high_freq_fft)
plt.figure(figsize=(6, 5))
plt.plot(time_vec, sig, label='Original signal')
plt.plot(time_vec, filtered_sig, linewidth=3, label='Filtered signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.legend(loc='best')
plt.show()
위의 코드는 다음과 같은 그래프를 결과로 출력합니다. 먼저 아래 그래프는 signal을 생성한 후 original signal을 표시한 것입니다.
이 original signal에 fft를 사용하여 power를 표시하면 다음과 같습니다.
모든 high frequencies를 제거하면 다음과 같은 그래프를 얻을 수 있습니다.
다음 예제를 실습해 보자. 포토샵 등에서 많이 사용하는 Gaussian image blur 효과를 사진에 적용하는 것입니다.
먼저 사진 이미지를 화면에 출력하는 코드를 실행 해 보겠습니다자.
import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt
# read image
img = plt.imread('./IMG_0820.png')
plt.figure()
plt.imshow(img)
plt.show()
IMG_0820.png 는 대전 동학사의 벗꽃 사진입니다. 위의 코드를 실행하면 다음과 같은 이미지가 나타납니다. 굉장히 크고 해상도가 좋은 사진 원본입니다.
이제 Gaussian convolution kernel 을 준비하고 FFT를 사용하여 적용해 봅니다.
전체 소스코드는 다음과 같습니다.
import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt
# read image
img = plt.imread('./IMG_0820.png')
plt.figure()
plt.imshow(img)
# First a 1-D Gaussian
t = np.linspace(-10, 10, 30)
bump = np.exp(-0.1*t**2)
bump /= np.trapz(bump) # normalize the integral to 1
# make a 2-D kernel out of it
kernel = bump[:, np.newaxis] * bump[np.newaxis, :]
# Padded fourier transform, with the same shape as the image
# We use :func:`scipy.signal.fftpack.fft2` to have a 2D FFT
kernel_ft = fftpack.fft2(kernel, shape=img.shape[:2], axes=(0, 1))
# convolve
img_ft = fftpack.fft2(img, axes=(0, 1))
# the 'newaxis' is to match to color direction
img2_ft = kernel_ft[:, :, np.newaxis] * img_ft
img2 = fftpack.ifft2(img2_ft, axes=(0, 1)).real
# clip values to range
img2 = np.clip(img2, 0, 1)
# plot output
plt.figure()
plt.imshow(img2)
plt.show()
