from scipy.signal import ellip, lfilter
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, lfilter_zi
def python_filtfilt(b, a, y, nfact=3, method='mirror'):
"""
Apply forward-backward filtering with initial conditions and bilateral padding.
Parameters:
b (ndarray): Numerator filter coefficients.
a (ndarray): Denominator filter coefficients.
y (ndarray): Input signal.
nfact (int, optional): Length of edge transients. Defaults to 3.
method ('pad' or 'mirror'): Padding method. Defaults to 'pad'.
Returns:
ndarray: Filtered signal.
"""
nfact=nfact*6
def create_mirrored_extension(yout, nfact):
"""
Create a mirrored extension for the given signal `yout`.
Parameters:
yout (ndarray): Original signal.
nfact (int): Length of the extension.
Returns:
ndarray: Mirrored extension.
"""
Npts = len(yout)
print(Npts)
# Left-side extension (linearly decreasing from yout[0] to yout[1])
left_ext = np.empty(nfact)
left_ext[0] = yout[0]
left_ext[1:] = 2 * yout[0] - yout[2:nfact + 1]
# Right-side extension (linearly increasing from yout[-1] to yout[-2])
right_ext = np.empty(nfact)
right_ext[-1] = yout[-1]
right_ext[:-1] = 2 * yout[-1] - yout[-2: -nfact - 1:-1]
# Concatenate the original signal with the extensions
ytemp = np.concatenate((left_ext, yout, right_ext))
print(ytemp.shape)
return ytemp
# Calculate initial conditions (zi) for lfilter
zi = lfilter_zi(b, a)
# Bilateral padding
Npts = len(y)
if method == 'pad':
# Pad with zeros on both ends
y_padded = np.concatenate((np.zeros(nfact), y, np.zeros(nfact)))
elif method == 'mirror':
# Mirror padding on both ends
y_padded = np.concatenate((
np.flip(2*y[1]-y[1:nfact+1], axis=0),
y,
np.flip(2*y[-1]-y[-nfact-2:-2], axis=0)
))
elif method=="linear":
y_padded = create_mirrored_extension(y, nfact)
else:
raise ValueError("Unsupported padding method. Choose 'pad' or 'mirror'.")
# Forward and Backward filtering
y_fwd, _ = lfilter(b, a, y_padded, zi=zi*y_padded[0])
y_rev, _ = lfilter(b, a, y_fwd[::-1], zi=zi*y_fwd[-1])
# Extract the central section of the filtered, reversed signal
y_out = y_rev[-(Npts+nfact):-nfact][::-1]
return y_out
data =np.load(r"2024-04-09 14_36_5_eegdata.npy")
signal =data[0,0,0,:]
# 现在 `filtered_signal` 是经过带通滤波处理后的信号数据
rp = 0.1 # 通带纹波（最大允许增益变化），一般取较小值如 0.1 或 1.0
rs = 80 # 阻带衰减（最小衰减量，以dB为单位），根据需要选择足够大的值以保证抑制效果
# 定义滤波器参数
fs = 250 # 采样率 (Hz)
lowcut = 6 # 通带下限 (Hz)
highcut = 48 # 通带上限 (Hz)
# 设计Butterworth带通滤波器
order =5 # 滤波器阶数，可以根据实际需求调整
nyq_freq = 0.5 * fs # Nyquist frequency
low = lowcut / nyq_freq
high = highcut / nyq_freq
b, a = butter(order, [low, high], btype='band')
y_filtered = python_filtfilt(b, a, signal, nfact=6, method='pad')
#将原始信号与滤波后的信号画图
#
plt.title("order:5 bandpassfilter:6-48")
plt.plot(signal, label='Original Signal')
plt.plot(y_filtered , label='y_filtered')
plt.legend()
plt.show()
