filiter.py

# TODO: PCG Processing
# - 使用环境变量/相对路径读取数据
# - 提供最小样例 data/PCG_ROW.csv
# FilterResponse.py
from  scipy import  signal
import matplotlib.pyplot as plt
import numpy as np

import pandas as pd
import warnings
import os

warnings.filterwarnings('ignore')

sr = 1000

start = 50000
stop = start+sr*5
data_path = os.getenv('PCG_ROW_PATH', 'data/PCG_ROW.csv')
data = pd.read_csv(data_path)
PCG = np.asarray(data)[start:stop]
def butterBandPassFilter(lowcut, highcut, samplerate, order):
	"生成巴特沃斯带通滤波器"
	semiSampleRate = samplerate * 0.5
	low = lowcut / semiSampleRate
	high = highcut / semiSampleRate
	b, a = signal.butter(order, [low, high], btype='bandpass')
	print("bandpass:", "b.shape:", b.shape, "a.shape:", a.shape, "order=", order)
	print("b=", b)
	print("a=", a)
	return b, a


def butterBandStopFilter(lowcut, highcut, samplerate, order):
	"生成巴特沃斯带阻滤波器"
	semiSampleRate = samplerate * 0.5
	low = lowcut / semiSampleRate
	high = highcut / semiSampleRate
	b, a = signal.butter(order, [low, high], btype='bandstop')
	print("bandstop:", "b.shape:", b.shape, "a.shape:", a.shape, "order=", order)
	print("b=", b)
	print("a=", a)
	return b, a

iSampleRate = 1000  # 采样频率,每秒2000个样本



plt.figure(figsize=(12, 5))
ax0 = plt.subplot(121)
for k in [2, 3, 4]:
	b, a = butterBandPassFilter(10, 50, samplerate=iSampleRate, order=k)
	w, h = signal.freqz(b, a, worN=2000)
	ax0.plot((iSampleRate * 0.5 / np.pi) * w, np.abs(h), label="order = %d" % k)

ax1 = plt.subplot(122)
for k in [2, 3, 4]:
	b, a = butterBandStopFilter(48, 52, samplerate=iSampleRate, order=k)
	w, h = signal.freqz(b, a, worN=2000)
	ax1.plot((iSampleRate * 0.5 / np.pi) * w, np.abs(h), label="order = %d" % k)

x = PCG
# 进行带通滤波
b, a = butterBandPassFilter(2,200, iSampleRate, order=4)
x = signal.lfilter(b, a, x)

# 进行带阻
b, a = butterBandStopFilter(49, 51, iSampleRate, order=2)
x = signal.lfilter(b, a, x)
plt.figure(figsize=(30,15))
plt.subplot(2,1,1)
plt.plot(PCG)
plt.subplot(2, 1, 2)
plt.plot(x)
plt.show()