Prolem с фильтром нижних частот в Python - PullRequest
Новая
       102

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

0 голосов
/ 28 августа 2018

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

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)

на который я получаю

this

с другой стороны, если я делаю то же самое, но с другим количеством очков, скажем, 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)

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

this

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

EDIT:

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

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))

на который я получаю this

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

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))

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

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

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

1 Ответ

0 голосов
/ 18 сентября 2018

Для тех, кого это может заинтересовать, на самом деле это «известная проблема» с масляным фильтром, используемым с выводом «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")

, что дает

this

...