코딩도상국

안녕하세요. 

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)

https://coding-yoon.tistory.com/23?category=830190

파이썬에는 장점 중 하나는 여러가지 라이브러리가 구현되어 있고, 손쉽게 이용할 수 있습니다.

위 글에 들어가보시면, Scipy 라이브러리에는 Low Pass Filter, High Pass Filter 등등 필터를 구현되어 있습니다. 

그렇다면 우리들 또한 라이브러리를 만들 수 있고, 배포할 수 있습니다.

이번 글은 Low Pass Filter, High Pass Filter, Band Pass Filter를 Scipy 라이브러리를 사용하지 않고 모듈화시키고 구현하겠습니다.

우선, 구현된 모습을 보겠습니다.

위 그림은 제가 Scipy 라이브러리를 사용하지 않았고, 제가 모듈화를 시켜서 돌린 것입니다.

파이썬에서 모듈화를 하는 방법은 굉장히 간단합니다.

예를 들어 filter_module의 이름을 가진 라이브러리라면, filter_module.py로 파이썬 파일을 만들어줍니다.

그리고 자신이 원하는 부분을 구현하면 됩니다.

단, 조건이 있다면 같은 프로젝트 안에 있어야 합니다. 아무리 같은 폴더 안에 있다 하더라도 같은 프로젝트로 공유하고 있지 않다면 인식을 할 수 없습니다. 

# filter_module.py
import numpy as np


# Low Pass Filter
def lpf_module(raw_data, time, f_cut_lpf, order, sr):
    """
    :param raw_data: Low Pass Filter 통과할 데이터
    :param time: 시간
    :param f_cut_lpf: 차단 주파수
    :param order: 차수
    :param sr: sample rate
    :param current_list: 진행도를 보여줄 list
    :return result: Low Pass Filter 데이터
    """
    sample_time = 1 / sr
    w_cut_lpf = 2 * np.pi * f_cut_lpf
    w_cut_period = 1 / w_cut_lpf
    raw_data_temp = raw_data.copy()
    previous_low_result = 0
    low_result = []
    low_list_result = []

    step = 0
    # low pass filter
    for k in range(order, 0, -1):
        for j in range(0, len(raw_data)):
            for i in range(0, len(time)):
                y = (w_cut_period * previous_low_result + sample_time * raw_data_temp[j][i]) / (
                        w_cut_period + sample_time)
                previous_low_result = y
                low_result.append(y)
            if k > 1:
                raw_data_temp[j] = low_result
                low_result = []
            if k == 1:
                low_list_result = np.append(low_list_result, low_result)
            low_result = []

    result = low_list_result.reshape(len(raw_data_temp), len(time))
    return result


# High Pass Filter
def hpf_module(raw_data, time, f_cut_hpf, order, sr):
    """
    :param raw_data: High Pass Filter 통과할 데이터
    :param time: 시간
    :param f_cut_hpf: 차단 주파수
    :param order: 차수
    :param sr: sample rate
    :param current_list: 진행도를 보여줄 list
    :return result: High Pass Filter 데이터
    """
    w_cut_hpf = 2 * np.pi * f_cut_hpf
    w_cut_period = 1 / w_cut_hpf
    sample_time = 1/sr
    """
    raw_data_temp = self.raw_data  # 주소를 복사
    raw_data_temp = self.raw_data.copy  # 값을 복사
    """
    raw_data_temp = raw_data.copy()  # raw data
    previous_high_result = 0  # previous result
    previous_raw_data = 0  # previous raw data
    high_result = []  # result
    high_list_result = []
    step = 0

    # high pass filter
    for k in range(order, 0, -1):
        for j in range(0, len(raw_data_temp)):
            for i in range(0, len(time)):
                y = w_cut_period / (w_cut_period + sample_time) * \
                    previous_high_result + w_cut_period / \
                    (w_cut_period + sample_time) * (
                            raw_data_temp[j][i] - previous_raw_data)
                previous_high_result = y
                previous_raw_data = raw_data_temp[j][i]
                high_result = np.append(high_result, y)
            if k > 1:
                raw_data_temp[j] = high_result
                high_result = []
            if k == 1:
                high_list_result = np.append(high_list_result, high_result)
            high_result = []

    # filter data reshape ( np.append 를 사용하면 일차배열로 쌓임 )
    result = high_list_result.reshape(len(raw_data_temp), len(time))
    return result


# Band Pass Filter
def bpf_module(raw_data, time, f_cut_hpf, f_cut_lpf, order, sr, current_list=None):
    """
        :param raw_data: Band Pass Filter 통과할 데이터
        :param time: 시간
        :param f_cut_lpf: 차단 주파수
        :param f_cut_hpf: 차단 주파수
        :param order: 차수
        :param sr: sample rate
        :param current_list: 진행도를 보여줄 list
        :return result: Band Pass Filter 데이터
        """
    sample_time = sr
    raw_data_temp = raw_data.copy()  # raw data
    w_cut_lpf = 2 * np.pi * f_cut_lpf
    w_cut_lpf_period = 1 / w_cut_lpf

    previous_low_result = 0  # previous raw data
    low_result = []  # result
    low_list_result = []
    previous_raw_data = 0  # previous raw data
    step = 0

    # low pass filter
    for k in range(order, 0, -1):
        for j in range(0, len(raw_data_temp)):
            for i in range(0, len(time)):
                y = (w_cut_lpf_period * previous_low_result + sample_time * raw_data_temp[j][i]) / (
                        w_cut_lpf_period + sample_time)
                previous_low_result = y
                low_result = np.append(low_result, y)
            if k > 1:
                raw_data_temp[j] = low_result
                low_result = []
            if k == 1:
                low_list_result = np.append(low_list_result, low_result)
            low_result = []

    raw_data_temp = low_list_result.reshape(len(raw_data_temp), len(time))
    previous_high_result = 0  # previous result
    w_cut_hpf = 2 * np.pi * f_cut_hpf
    w_cut_period = 1 / w_cut_hpf
    high_result = []  # result
    high_list_result = []
    step = 0
    # high pass filter
    for k in range(order, 0, -1):
        for j in range(0, len(raw_data_temp)):
            for i in range(0, len(time)):
                y = w_cut_period / (w_cut_period + sample_time) * \
                    previous_high_result + w_cut_period / \
                    (w_cut_period + sample_time) * (
                            raw_data_temp[j][i] - previous_raw_data)
                previous_high_result = y
                previous_raw_data = raw_data_temp[j][i]
                high_result = np.append(high_result, y)
            if k > 1:
                raw_data_temp[j] = high_result
                high_result = []
            if k == 1:
                high_list_result = np.append(high_list_result, high_result)
            high_result = []
    # filter data reshape ( np.append 를 사용하면 일차배열로 쌓임 )
    result = high_list_result.reshape(len(raw_data_temp), len(time))
    return result

이렇게 Low Pass Filter, High Pass Filter, Band Pass Filter를 구현된 저만의 라이브러리를 만들었습니다. 

그렇다면 어떻게 사용할까요? 사용방법도 굉장히 간단합니다. 평소 우리들이 다운받아 사용하는 라이브러리처럼

import 해서 불러오면 됩니다. from에는 파일명이 들어갑니다. 

from filter_module import *
"""
from filter_module import lpf_module
from filter_module import hpf_module
from filter_module import bpf_module
"""

raw_data = "데이터가 담겨 있는 이차원 배열형식"
time = "데이터의 시간영역"
cut_off_frequency = "차단주파수"
order = "차수"
sr = "sample rate"

# 저역통과필터
lpf_result = lpf_module(raw_data, time, cut_off_frequency, order, sr)
# 고역통과필터
hpf_result = hpf_module(raw_data, time, cut_off_frequency, order, sr)

이런식으로 모듈화를 시켜 프로그래밍을 진행하게 되면, 나중에 팀프로젝트할 때도 상당히 수월하며, 몇 달이 지나 수정

하고 싶다고 해도 비교적 쉽게 할 수 있습니다. 그리고 모듈화의 큰 장점은 재사용성입니다. 

나중에 필터기능이 필요할 때, 위 코딩처럼 주석과 모듈화를 시켜놓으면 filter_module.py만 불러오기만 하면 됩니다. 

처음엔 복잡할지라도 점점 코딩이 이뻐지기 시작하면서 기분도 좋아집니다@@

안녕하세요. 글은 계속 쓰는데 적는데 마음에 안들어서 전부 비공개로 해놨는데...

네 그렇다구요.

이번 글은 파이썬에서 Scipy를 통해서 음성신호를 필터링해볼려고 합니다. 

음성신호를 필터링했으면 matplotlib 라이브러리로 파형 한 번 확인해보고,

샘플링 된 음성파일을 저장해 들어볼려고 합니다.

https://www.scipy.org/

위 주소는 Scipy, 사이파이라고 과학기술계산을 위한 Python 라이브러리입니다.  

우선 그럼 설치부터 해볼까요.

cmd 창에 아래처럼 적어줍시다.

pip install scipy

만약 Fatal error in launcher : 블라블라 라고 뜨면

python -m pip install scipy

이렇게 적어줍니다. 아마 아나콘다 다운로드를 하면 자동적으로 다운되는 걸로 알고 있습니다.

다운로드가 되었다면 import 해봅시다.

import가 잘되네요. 

음성 파일 경로, 이름 한 번씩 확인해볼까요?

scipy.io => io : input, output

경로는 data에 들어가 있고, 이름은 test네요.

음성 파일 불러온 후에 파형 한 번 확인해봅시다.

wavefile.read는 scipy.io에 들어있습니다. 

x축으로는 시간, y축으로는 진폭을 나타냅니다. 

혹시 모르니 sr, x 값 찍어볼까요?

sr = sample rate는 48000이고

x = wav의 데이터들이 list로 들어가 있습니다. 

파형을 보니 음성파일 시간은 3초짜리 영상이네요. 

(오늘은 scipy에 대해 알아보는 시간이니 numpy, matpltlib는 생략하도록 하겠습니다. )

그럼 먼저 LPF(Low Pass Filter : 저역통과필터) 한 번 돌려볼까요?

오늘의 핵심 포인트!!

signal.firwin 이 함수 하나로 hpf, lpf, bpf 다 구현할 수 있습니다. 함수 원형 한 번 확인해볼까요.

https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.firwin.html

signal.firwin(101, cutoff=500, fs=sr, pass_zero='lowpass')

함수 중요한 부분만 적어본다면,

signal.firwin(차수, 차단주파수, 필터)

차수 : 영어로 order, 필터의 길이입니다. 어느정도 값을 크게 줘야합니다. 

차단주파수 : 어디서부터 주파수를 자를 것인지 결정합니다. 

필터 : 말보단 사진으로 보여드리겠습니다.

오... bandpass, lowpass, highpass 뿐만아니라 bandstop도 되는군요. 

만약 쓰고 싶은 필터가 있으면 pass_zero = 'bandstop' 이렇게 하시면 됩니다.

제가 적용한 함수를 보면

차수는 101, 차단주파수는 500hz, 필터는 low pass filter를 사용합니다. 

주의사항으로 차수는 must be odd if a passband includes the Nyquist frequency라고

반드시,  나이키스트 주파수 때문에 홀수를 사용해야 한다고 합니다.

파형 보겠습니다.

아래 주파수는 잘린 것을 확인할 수 있습니다.

그럼 이제 low pass 필터링된 음성들어보겠습니다.

잘 들어가 있네요?

혹여나 경로가 이상하다고 생각하실 수 있는데

저 같은 경우는 주피터 노트북을 실행할 때

cd C:\Users\Mangnani\Desktop\train

변경을 해줬기 때문에 기본 주소가 train으로 잡혀있습니다. 

요렇게요. 그럼 한 번 들어보겠습니다.

네. 필터링이 아주 잘된 것 같습니다. 하하..잘들리시나요?

그럼 High Pass Filter 가보겠습니다.

이번에는 한 번에 적었습니다. 한 번 들어볼까요?

네.. 이번에도 잘된 것 같습니다.

만약 필터링을 해보고 싶으시다면 위와 같이 하시면 될 것 같습니다.

좀 더 자세하게 깊게 하고 싶으시다면 차수, 차단주파수 등등 위에 있는 레퍼런스에 따라 

실험해보시면 될 것 같습니다.