|
from scipy.signal import butter, lfilter
|
|
import torch
|
|
from scipy import signal
|
|
import librosa
|
|
import numpy as np
|
|
|
|
from scipy.signal import sosfiltfilt
|
|
from scipy.signal import butter, cheby1, cheby2, ellip, bessel
|
|
from scipy.signal import resample_poly
|
|
|
|
|
|
def align_length(x=None, y=None, Lx=None):
|
|
"""align the length of y to that of x
|
|
|
|
Args:
|
|
x (np.array): reference signal
|
|
y (np.array): the signal needs to be length aligned
|
|
|
|
Return:
|
|
yy (np.array): signal with the same length as x
|
|
"""
|
|
assert y is not None
|
|
|
|
if Lx is None:
|
|
Lx = len(x)
|
|
Ly = len(y)
|
|
|
|
if Lx == Ly:
|
|
return y
|
|
elif Lx > Ly:
|
|
|
|
return np.pad(y, (0, Lx - Ly), mode="constant")
|
|
else:
|
|
|
|
return y[:Lx]
|
|
|
|
|
|
def bandpass_filter(x, lowcut, highcut, fs, order, ftype):
|
|
"""process input signal x using bandpass filter
|
|
|
|
Args:
|
|
x (np.array): input signal
|
|
lowcut (float): low cutoff frequency
|
|
highcut (float): high cutoff frequency
|
|
order (int): the order of filter
|
|
ftype (string): type of filter
|
|
['butter', 'cheby1', 'cheby2', 'ellip', 'bessel']
|
|
|
|
Return:
|
|
y (np.array): filtered signal
|
|
"""
|
|
nyq = 0.5 * fs
|
|
lo = lowcut / nyq
|
|
hi = highcut / nyq
|
|
|
|
if ftype == "butter":
|
|
|
|
sos = butter(order, [lo, hi], btype="band", output="sos")
|
|
elif ftype == "cheby1":
|
|
sos = cheby1(order, 0.1, [lo, hi], btype="band", output="sos")
|
|
elif ftype == "cheby2":
|
|
sos = cheby2(order, 60, [lo, hi], btype="band", output="sos")
|
|
elif ftype == "ellip":
|
|
sos = ellip(order, 0.1, 60, [lo, hi], btype="band", output="sos")
|
|
elif ftype == "bessel":
|
|
sos = bessel(order, [lo, hi], btype="band", output="sos")
|
|
else:
|
|
raise Exception(f"The bandpass filter {ftype} is not supported!")
|
|
|
|
|
|
y = sosfiltfilt(sos, x)
|
|
|
|
if len(y) != len(x):
|
|
y = align_length(x, y)
|
|
return y
|
|
|
|
|
|
def lowpass_filter(x, highcut, fs, order, ftype):
|
|
"""process input signal x using lowpass filter
|
|
|
|
Args:
|
|
x (np.array): input signal
|
|
highcut (float): high cutoff frequency
|
|
order (int): the order of filter
|
|
ftype (string): type of filter
|
|
['butter', 'cheby1', 'cheby2', 'ellip', 'bessel']
|
|
|
|
Return:
|
|
y (np.array): filtered signal
|
|
"""
|
|
nyq = 0.5 * fs
|
|
hi = highcut / nyq
|
|
|
|
if ftype == "butter":
|
|
sos = butter(order, hi, btype="low", output="sos")
|
|
elif ftype == "cheby1":
|
|
sos = cheby1(order, 0.1, hi, btype="low", output="sos")
|
|
elif ftype == "cheby2":
|
|
sos = cheby2(order, 60, hi, btype="low", output="sos")
|
|
elif ftype == "ellip":
|
|
sos = ellip(order, 0.1, 60, hi, btype="low", output="sos")
|
|
elif ftype == "bessel":
|
|
sos = bessel(order, hi, btype="low", output="sos")
|
|
else:
|
|
raise Exception(f"The lowpass filter {ftype} is not supported!")
|
|
|
|
y = sosfiltfilt(sos, x)
|
|
|
|
if len(y) != len(x):
|
|
y = align_length(x, y)
|
|
|
|
y_len = len(y)
|
|
|
|
y = stft_hard_lowpass(y, hi, fs_ori=fs)
|
|
|
|
y = sosfiltfilt(sos, y)
|
|
|
|
if len(y) != y_len:
|
|
y = align_length(y=y, Lx=y_len)
|
|
|
|
return y
|
|
|
|
|
|
def stft_hard_lowpass(data, lowpass_ratio, fs_ori=44100):
|
|
fs_down = int(lowpass_ratio * fs_ori)
|
|
|
|
y = resample_poly(data, fs_down, fs_ori)
|
|
|
|
|
|
y = resample_poly(y, fs_ori, fs_down)
|
|
|
|
if len(y) != len(data):
|
|
y = align_length(data, y)
|
|
return y
|
|
|
|
|
|
def limit(integer, high, low):
|
|
if integer > high:
|
|
return high
|
|
elif integer < low:
|
|
return low
|
|
else:
|
|
return int(integer)
|
|
|
|
|
|
def lowpass(data, highcut, fs, order=5, _type="butter"):
|
|
"""
|
|
:param data: np.float32 type 1d time numpy array, (samples,) , can not be (samples, 1) !!!!!!!!!!!!
|
|
:param highcut: cutoff frequency
|
|
:param fs: sample rate of the original data
|
|
:param order: order of the filter
|
|
:return: filtered data, (samples,)
|
|
"""
|
|
|
|
if len(list(data.shape)) != 1:
|
|
raise ValueError(
|
|
"Error (chebyshev_lowpass_filter): Data "
|
|
+ str(data.shape)
|
|
+ " should be type 1d time array, (samples,) , can not be (samples, 1)"
|
|
)
|
|
|
|
if _type in "butter":
|
|
order = limit(order, high=10, low=2)
|
|
return lowpass_filter(
|
|
x=data, highcut=int(highcut), fs=fs, order=order, ftype="butter"
|
|
)
|
|
elif _type in "cheby1":
|
|
order = limit(order, high=10, low=2)
|
|
return lowpass_filter(
|
|
x=data, highcut=int(highcut), fs=fs, order=order, ftype="cheby1"
|
|
)
|
|
elif _type in "ellip":
|
|
order = limit(order, high=10, low=2)
|
|
return lowpass_filter(
|
|
x=data, highcut=int(highcut), fs=fs, order=order, ftype="ellip"
|
|
)
|
|
elif _type in "bessel":
|
|
order = limit(order, high=10, low=2)
|
|
return lowpass_filter(
|
|
x=data, highcut=int(highcut), fs=fs, order=order, ftype="bessel"
|
|
)
|
|
|
|
|
|
|
|
|
|
else:
|
|
raise ValueError("Error: Unexpected filter type " + _type)
|
|
|
|
|
|
def bandpass(data, lowcut, highcut, fs, order=5, _type="butter"):
|
|
"""
|
|
:param data: np.float32 type 1d time numpy array, (samples,) , can not be (samples, 1) !!!!!!!!!!!!
|
|
:param lowcut: low cutoff frequency
|
|
:param highcut: high cutoff frequency
|
|
:param fs: sample rate of the original data
|
|
:param order: order of the filter
|
|
:param _type: type of filter
|
|
:return: filtered data, (samples,)
|
|
"""
|
|
if len(list(data.shape)) != 1:
|
|
raise ValueError(
|
|
"Error (chebyshev_lowpass_filter): Data "
|
|
+ str(data.shape)
|
|
+ " should be type 1d time array, (samples,) , can not be (samples, 1)"
|
|
)
|
|
if _type in "butter":
|
|
order = limit(order, high=10, low=2)
|
|
return bandpass_filter(
|
|
x=data,
|
|
lowcut=int(lowcut),
|
|
highcut=int(highcut),
|
|
fs=fs,
|
|
order=order,
|
|
ftype="butter",
|
|
)
|
|
elif _type in "cheby1":
|
|
order = limit(order, high=10, low=2)
|
|
return bandpass_filter(
|
|
x=data,
|
|
lowcut=int(lowcut),
|
|
highcut=int(highcut),
|
|
fs=fs,
|
|
order=order,
|
|
ftype="cheby1",
|
|
)
|
|
|
|
|
|
elif _type in "ellip":
|
|
order = limit(order, high=10, low=2)
|
|
return bandpass_filter(
|
|
x=data,
|
|
lowcut=int(lowcut),
|
|
highcut=int(highcut),
|
|
fs=fs,
|
|
order=order,
|
|
ftype="ellip",
|
|
)
|
|
elif _type in "bessel":
|
|
order = limit(order, high=10, low=2)
|
|
return bandpass_filter(
|
|
x=data,
|
|
lowcut=int(lowcut),
|
|
highcut=int(highcut),
|
|
fs=fs,
|
|
order=order,
|
|
ftype="bessel",
|
|
)
|
|
else:
|
|
raise ValueError("Error: Unexpected filter type " + _type)
|
|
|