Prolem с фильтром нижних частот в Python

Question

Я хотел отфильтровать (низкочастотный) сигнал, который у меня есть, и когда он не работал, я начал исследовать, почему это не будет. Я сделал несколько тестов, и я несколько удивлен поведением фильтра Баттерворта. я определил это как в этом посте

def apply_filter(data, cutoff, fs, order=6, filter_type="low", analog=False):
    nyq = 0.5 * fs
    normalized_cutoff = cutoff / nyq
    b,a = butter(order, normalized_cutoff, btype=filter_type, analog=analog, output="ba")
    they = lfilter(b, a, data)
    return(they)

если я возьму образец длиной 1000 элементов, вот так

x = np.linspace(0, 2*np.pi, 1000)
y = np.sin(x) + 0.3* np.sin(10*x)
sampling_frequency = 1/ (x[-1] * 1e-3)
sampling_frequency
>> 159.15494309189532
# because i have 1000 thousand points for a "time" going up to 2 pi

plt.plot(x, y, x, apply_filter(y, cutoff=1, fs= sampling_frequency)

на что я получаю

с другой стороны, если я делаю то же самое, но с другим количеством очков, скажем, 10000, я получаю неправильный результат, и я не совсем понимаю, почему:

x = np.linspace(0, 2*np.pi, 10000)
y = np.sin(x) + 0.3* np.sin(10*x)
sampling_frequency = 1/ (x[-1] * 1e-4)
sampling_frequency
>> 1591.5494309189535
# because i have 10000 thousand points for a "time" going up to 2 pi

plt.plot(x, y, x, apply_filter(y, cutoff=1, fs= sampling_frequency)

и на этот раз я получаю

что явно не так. Может кто-нибудь объяснить, почему это произошло? кажется, что все работает нормально на 1000 очков или меньше...

РЕДАКТИРОВАТЬ:

Я построил частотную характеристику фильтра, и проблема появляется на этих графиках, хотя я тоже не знаю, почему это происходит.

sampling rate 
>> 159.1549430918953
b, a = butter(6, 1/(sampling_rate/2))
w, h = freqz(b, a, 8000)
plt.subplot(2,1,1)
plt.xlim(0, 15)
plt.plot(0.5*sampling_rate*w/np.pi, np.abs(h))

на что я получаю

тогда как, если я сделаю

sampling_frequency *= 10
sampling_frequency
>> 1591.5494309189535
b, a = butter(6, 1/(sampling_rate/2))
w, h = freqz(b, a, 8000)
plt.subplot(2,1,1)
plt.xlim(0, 15)
plt.plot(0.5*sampling_rate*w/np.pi, np.abs(h))

тогда я получаю

Я чувствую, что у функции Баттерворта есть какие-то проблемы с большим количеством очков по какой-то причине?

Спасибо за вашу помощь!

Answer 1

Для тех, кого это может заинтересовать, на самом деле это "известная проблема" с масляным фильтром, используемым с выводом "ba". Вместо этого вы должны использовать вывод "zpk", см. Эту ссылку.

Вы можете использовать вывод "zpk" довольно простым способом, очень похожим на то, что вы делаете с выводом "ba". Похоже, это работает на 1 миллион очков, и нет никаких причин, чтобы не работать дальше.

вот основной пример:

point_number=1000000
# our "data"
x = np.linspace(0, 2*np.pi, point_number)
y = sin(x) + 0.3* sin(10*x)

# sampling frequency would be 1/ sampling_time
sampling_frequency = point_number/(2*np.pi)
# hence the nyquist frequency
nyq = sampling_frequency/2
# desired cutoff frequency, in Hertz
cutoff = 1
# normalized for the function butter
normalized_cutoff = cutoff/nyq

z,p,k = butter(6, normalized_cutoff, output="zpk")
lesos = zpk2sos(z, p, k)
# filtered data
y_filtered_sos = sosfilt(lesos, y)

# plot part
plt.plot(x, y_filtered_sos, label="filtered data")
plt.title("filtered data, with zpk")
plt.plot(x,y, label="data")
plt.legend()
plt.title("filtered data, with zpk")

который дает