์๋ ํ์ธ์.
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)
728x90
๋ฐ์ํ