코딩도상국

안녕하세요. 

Low Pass Filter, High Pass Filter, Band Pass Filter를 이용한 신호 처리와 FFT 변환으로 주파수 영역에서 필터가 제대로 동작했는지 확인하는 것까지 파이썬으로 구현합니다. 

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
import scipy.io
import os


mat_file = scipy.io.loadmat('signal1.mat')
(file_path, file_id) = os.path.split('signal1.mat')  # file path, file name

fs = 1024  # sample rate
order = 10  # order
cut_off_freq = 150  # cut off frequency

# raw signal
# t = np.linspace(0, 1, fs, False)  # 1 second
# sig = np.sin(2*np.pi*100*t) + np.sin(2*np.pi*50*t)  # signal
sig = mat_file[file_id[:-4]][0]

freq = np.fft.fftfreq(len(sig), 1/1024)

# filtered signal
sos = signal.butter(order, [cut_off_freq], 'low', fs=fs, output='sos')  # low pass filter
filtered = signal.sosfilt(sos, sig)

# raw signal fft
raw_fft = np.fft.fft(sig) / len(sig)
raw_fft_abs = abs(raw_fft)

# filter signal fft
filtered_fft = np.fft.fft(filtered) / len(filtered)
filtered_fft_abs = abs(filtered_fft)

# plot
fig, ((ax00, ax01), (ax10, ax11)) = plt.subplots(2, 2)

# raw signal plot : 0 row 0 column
ax00.plot(t, sig)
ax00.set_title('Raw Data Time Domain')
ax00.set_xlabel('Time [seconds]')
ax00.set_ylabel('Amplitude')

# filtered signal plot : 1 row 0 column
ax10.plot(t, filtered)
ax10.set_title('Filtered Data Time Domain')
ax10.set_xlabel('Time [seconds]')
ax10.set_ylabel('Amplitude')

# raw signal fft plot : 0 row 1 column
ax01.stem(freq, raw_fft_abs, use_line_collection=True)
ax01.set_title('Raw Data Frequency Domain')
ax01.set_xlabel('Frequency [HZ]')
ax01.set_ylabel('Amplitude')

# filtered signal fft plot : 1 row column
ax11.stem(freq,filtered_fft_abs, use_line_collection=True)
ax11.set_title('Filtered Data Frequency Domain')
ax11.set_xlabel('Frequency [HZ]')
ax11.set_ylabel('Amplitude')

# plot
plt.tight_layout()
plt.show()

1. Low Pass Filter

sos = signal.butter(order, [cut_off_freq], 'low', fs=fs, output='sos')
filtered = signal.sosfilt(sos, sig)

2. High Pass Filter

sos = signal.butter(order, [cut_off_freq], 'high', fs=fs, output='sos')
filtered = signal.sosfilt(sos, sig)

3. Band Pass Filter

sos = signal.butter(order, [150, 200], 'band', fs=fs, output='sos')
filtered = signal.sosfilt(sos, sig)